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Chapter 1 
Introduction 



Mathematical modelling is that abstruse subject which forms the connecting tissue 
between the problems of the real world which we wish to solve, and the quantitative 
analysis which we undertake to do so. Almost any problem which requires a quanti- 
tative answer, whether it be in industry, medicine, economics, biology or geophysics 
(for example) involves the formulation of the problem as a mathematical model, and 
it is this formulation, and the techniques which one uses to solve the model, which 
form the subject of these notes. 

There are, perhaps, three kinds of model: statistical, discrete and continuous. 
We can illustrate their difference with a simple (and topical) example. Suppose fox 
hunting in England were to be banned. One of the arguments against this is that 
hunting provides a control on the number of foxes. As with many such assertions, 
this is one which might be true, or might not. It is a quantitative assertion, and the 
only way to find out whether it is true is to examine it scientifically. We do this by 
proposing a model, and examining its validity 1 . The statistical method examines the 
probability that one of two alternative hypotheses is true, by using actual measured 
data. For example, if hunting were to be banned, one might examine fox population 
numbers during the five years preceding and following such a ban. As with all models, 
it is in the interpretation of the results that the skill lies. This is particularly true of 
statistics, where the aim is to eliminate possibilities, rather than propose constituent 
mechanisms. 

Statistical models have to deal with the issue of predictability. Suppose fox num- 
bers double following a ban; then it would seem likely that there is a connection. But 
there are other factors which are involved in determining fox populations, and it is 
always arguable whether one particular factor has a deciding control. 

Statistical models are diagnostic: they try to intepret process from measured data. 
Discrete and continuous models are examples of prognostic models: they propose a 
descriptive model of a phenomenon, and then predict what will happen in the future. 
Such models require validation. This consists initially of matching observation to 
theory, and often will suggest experiments which can be done to confirm the theory. 

A discrete model will propose a difference equation for the variables of interest. 

J This is in fact the scientific method: the models that we propose are no more than hypotheses, 
and science should never present itself as purveying absolute truth. 
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In the case of foxes, this might be an estimate for the fox population density u n at a 
particular time during the n-th year, and a simple such discrete model is the logistic 
model with harvesting: 

u n+1 = ru n (l - j^J - hu n , (1.1) 

where r is the specific reproductive rate, and the term hu n represents harvesting (via 
hunting); the coefficient h represents the effort involved, which one might suppose 
would be proportional to the number of hunts. This model contains the simple idea 
that excess populations become limited by competition for resources (the nonlinear 
term involving K implies decreasing growth rate at larger values of u n ). A discrete 
model such as (1.1) might be appropriate for fox populations, which have an annual 
rut, so that the reproductive cycle is repeated annually. 

Suppose, for the sake of argument, that (1.1) is a reasonable model. Note that 
r > 1 is a pure number, and so also is h, while K has the same units as u n . In the 
absence of hunting, there is a steady population K(r — l)/r, and this is reduced if 

. So the efficacy of hunting in this model depends on 

r — 1J 

this particular quantity; hunting is effective if h/(r — 1) is significant. This provides 
a predicted outcome, provided the parameters h and r can be reasonably assessed. 

Continuous models are used in the same way, but describe processes using dif- 
ferential equations. It is this kind of model which forms the focus of these notes. 
The rationale for continuous models is the continuum hypothesis, which states that 
the actual behaviour of a discrete variable, such as a population, can be accurately 
represented by the evolution of a continuous (and usually differentiable) variable. 
The basis for this assumption is the 'fine-grained' nature of the variable: in a large 
population, a change by one individual is a small relative change, and can be viewed 
as being a finite difference approximation to the infinitesimal changes of the calculus. 

A continuous version of (1.1) might be the ordinary differential equation 

_ = /ra ( 1 __j_„„, (1 . 2) 

whose behaviour can be analysed in a similar way to its discrete counterpart. 

There are two particular points of view which we can bring to bear on the math- 
ematical models which describe the phenomena which concern us: these are the dy- 
namical systems approach, or equivalently the bifurcation theory approach; and the 
perturbation theory approach. Each has its place in different contexts, and sometimes 
they overlap. 

The bifurcation theory approach most usually (but not always) is brought to bear 
on problems which have some kind of complicated time-dependent behaviour. The 
idea is that we seek to understand the observations through the understanding of 
a number of simpler problems, which arise successively through bifurcations in the 
mathematical model, as some critical parameter is changed. A classic example of this 
approach is in the study of the origin of chaos in the Lorenz equations, or the onset 
of complicated forms of thermal convection in fluids. 
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In its simplest form (e.g., in weakly nonlinear stability theory) the perturbative 
approach is similar in method to the bifurcational one; however, the ethos is rather 
different. Rather than try and approach the desired solution behaviour through a se- 
quence of simpler behaviours, we try and break down the solution by making approxi- 
mations, which (with luck) are in fact realistic. In real problems, such approximations 
are readily available, and part of the art of the applied mathematician is having the 
facility of being able to judge how to make the right approximations. In these notes, 
we follow the perturbative approach. It has the disadvantage of being harder, but it 
is able to get closer to a description of how realistic systems may actually behave. 



1.1 Conservation and constitutive laws 

The basic building blocks of continuous mathematical models are conservation laws. 
The continuum assumption adopts the view that the physical medium of concern 
may be considered continuous, whether it be a porous medium (for example, sand 
on a beach) or a fluid flow. The continuum hypothesis works whenever the length 
or time scales of interest are (much) larger than the corresponding microscale. For 
example, the formation of dunes in a desert (length scale hundreds of metres) can 
be modelled as a continuous process, since the microscale (sand grain size) is much 
smaller. Equally, the modelling of large animal populations or of snow avalanches 
treats the corresponding media as continuous. 

Conservation laws arise as mathematical equations which represent the idea that 
certain quantities are conserved — for example, mass, momentum (via Newton's law) 
and energy. More generally, a conservation law refers to an equation which relates 
the increase or decrease of a quantity to terms representing supply or destruction. 

In a continuous medium, the typical form of a conservation law is as follows: 

g + V.f = * (1.8) 

In this equation, <p is the quantity being 'conserved' (expressed as amount per unit 
volume of medium, i.e. as a density; f is the 'flux', representing transport of cf) within 
the medium, and S represents source (S > 0) or sink (S < 0) terms. Derivation of 
the point form (1.3) follows from the integral statement 



d 
d~t 



f(j>dV = -f f.ndS+ [ SdV, (1.4) 

JV JdV JV 



after application of the divergence theorem (which requires f to be continuously dif- 
ferentiate), and by then equating integrands, on the basis that they are continuous 
and V is arbitrary. Derivation of (1.3) thus requires cf) and f to be continuously 
differentiable, and S to be continuous. 

Two basic types of transport are advection (the medium moves at velocity u, so 
there is an advective flux </>u) and diffusion, or other gradient-driven transport (such 
as chemotaxis). One can thus write 

f = 0u + J, (1.5) 
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where J might represent diffusive transport, for example. The very simplest conser- 
vation law is that of conservation of mass, where the conserved quantity is the density 
p, and the mass flux is entirely due to advection: 

^ + V.(pu) = 0. (1.6) 

Invariably, conservation laws give more terms than equations. In (1.5), for exam- 
ple, we have one scalar equation for </>, but other quantities J and S are present as 
well, and equations for these must be provided. Typically, these take the form of con- 
stitutive laws, and are based squarely on experimental measurement. For example, 
diffusive transport is represented by the assumption 

J = -DV(f>, (1.7) 

where D is a diffusion coefficient. In the heat equation, this is known as Fourier's 
law, and the heat equation itself takes the familiar form 

j t {pc P T) + V.[pc p Tu] = V.[kVT] + Q, (1.8) 
where Q represents any internal heat source or sink. 



1.2 The law of mass action 

Apart from conservation laws, which typically involve transport terms in the form of 
divergences, the other type of term which frequently arises in models is an algebraic 
'reaction' term. The source term Q in (1.8) is a simple example of this. Slightly more 
complicated is the Newtonian cooling law 

^ = h(T -T), (1.9) 

which represents cooling (or warming) of a material of temperature T in surroundings 
of ambient temperature T . Here, hT represents heat gain from the surroundings, 
and hT is the corresponding heat loss. 

A particularly common and interesting way in which source terms arise is in the 
nonlinear interaction between different constituent populations. These may represent 
individuals of a plant or animal species, or they may represent different types of 
cells in the body, or concentrations of different reacting chemical species. The basis 
for describing such interactions is called the law of mass action, and is most simply 
understood in the reaction between two substances, which we denote as A and B. 
The reaction is written as 

A + B^P, (1.10) 

indicating that a product P is formed at a rate r, which may typically have dimensions 
mol 1 _1 s _1 ; 'mol' indicates moles, and '1' indicates litres. The law of mass action says 
that 

roc AB, (1.11) 
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where A and B indicate the concentrations of the respective reactants. This law 
represents the idea that the rate of a chemical reaction is proportional to the rate at 
which the reacting molecules run into each other, and in a homogeneous, well-mixed 
environment, this is simply proportional to the size of each population. Insofar as the 
same basic assumption applies to any population, the law of mass action will apply 
there too. 



1.3 The S-I-R model 

A simple illustration of the law of mass action is the SIR model for epidemics in 
a population of fixed size. The population is divided into three classes, susceptible 
(those who have yet to catch the disease), infected, and removed (or recovered). 
The removed class may represent those who have recovered from the disease and are 
no longer susceptible (as, for example, in a flu epidemic) or those who have died 
from the disease. Three equations which describe the interactions of these three 
sub-populations are 

S = -kSI, 
I = kSI-rl, 

R = rl. (1.12) 

The nonlinear interaction term kSI represents the rate of infection, and assumes 
free dispersal of the infected class amongst the general population. Also of note is 
the first order decay term —rl, which describes 'first order' kinetics, and is indicative 
of the law of mass action applied to a single population. However, in the case of an 
epidemic, this term is of dubious provenance, as it implies that there is an exponential 
distribution of infection times, whereas in reality diseases tend to last for a reasonably 
fixed period. Consideration of this leads us to consider an age-structured model. 

Let us suppose that i(t,a) represents the density (with respect to a) of infected 
individuals at time t who have been infected for a time a (i. e., their 'age' of infection 
is a). We suppose the rate of removal r(a) is a function of age, so that 

di 

- = -r(a)i, (1.13) 
and the evolution of age with time is simply 

These are the characteristic equations for the age-structured evolution equation 

at + aa = - r(o)1 ' (L15) 



The total infected class is 



f oo 

I 



= ida, (1.16) 
Jo 
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which on integrating (1.15) implies 

dl r°° 

Evidently i\ a=0 is the rate of addition to the infected class, thus we take 

i\ a=0 = kSI, (1.18) 

and this provides the initial condition for i on a = 0. 

If r is constant we recover (1.12). For r = r(a), we solve (1.15) using the method 
of characteristics, i.e., we solve the ordinary differential equations (1.13) and (1.14) 
subject to the initial conditions on a = which can be parameterised as 

t = v , i = io(v) = kS(ri)I(ri), a = 0, (1.19) 

valid for rj > 0. (The other part of the initial condition is applied at t = 0, a > 0, but 
this part of the initial condition only contributes to a transient part of the solution 
for a finite time.) 

The solution of (1.13) and (1.14) is then found to be 



i — io(t — a) exp 



ra 

f r(s) 
Jo 



ds 



(1.20) 



valid for t > a. For t < a, the solution depends on the initial condition at t = 0, and 
if, for example, i — initially, then i = for t < a. 

Now let us consider the particular case where there is a fixed period of infection, 
r. This means essentially that r = for a < r, but becomes rapidly infinite 2 as a 
increases through r. As a consequence, we have that 



and thus 

i — io(t — a), a < r, 
% = 0, a > t. 



(1.22) 



Finally we have that (for t > r) 

1=1* i ( s )ds. (1.23) 

Jt-T 

In view of the definition of io, we can deduce from this that 

I = -S(t) + S(t-T), (1.24) 
and therefore that S satisfies the differential delay equation 

S = kS[S-S T ], (1.25) 

where S T = S(t — r). 

2 r is actually then what is called a generalised function, or distribution. 
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1.4 Saturation and cooperativity 



The simple quadratic interaction terms of the law of mass action become more compli- 
cated either when there are more reactions, or when the reaction rates are themselves 
modified by external factors. An example of the latter sort of term arises in the simple 
model for spruce budworm infestations given by Murray (2002). The model describes 
the evolution of spruce budworm moth density, denoted by B, as 



The spruce budworm larvae develop in pine forests in Northern Canada, where 
they feed on the foliage. The model describes the growth of the larval density as a 
term proportional to population, tbB (representing the year to year growth of larval 
populations through the successive stages of pupation, moth development and egg 
formation). The specific growth rate tb is modified by a term which represents the 
effects of crowding on the population. Thus we have the logistic growth model 



The quantity K B represents a saturation limit beyond which B cannot grow, and this 
represents one way in which the effects of saturation may be modelled. It is not the 
only way, and not necessarily the most realistic. 

The second term on the right hand side represents predation of the larvae by 
birds, for whom the larvae form part (but certainly not all) of the staple diet. The 
form of this term is of some interest. It is apparently independent of the size of the 
bird population. This implicitly assumes that there are so many birds that there is 
always one available to eat available larvae. Alternatively, it assumes a constant bird 
population, whose size is not controlled by larval density, but by other factors such 
as climate. 

At high values of B, the predation rate saturates. Presumably this represents 

the fact that the appetite of birds and other predators is finite. Saturation is often 

B 

modelled by a logistic term oc — ^ , which is a linear at small B, corresponding to 

first order kinetics. When the dependence is quadratic, as in (1.26), the rate is said 
to display the Allee effect. There are (at least) two motivations for such a quadratic 
dependence. The first is that when B is very small, one can imagine birds finding other 
prey more easily, and making less effort to locate larvae. A quadratic dependence 
would indicate that the effort or time spent searching for larvae is proportional to their 
density. The second is that at higher densities, the availability of larvae may attract 
birds. If bird density is proportional to larval density, then a quadratic dependence 
also ensues. 

1.4.1 Cooperativity 

Cooperativity refers to a biochemical reaction process in which an enzyme reacts with 
a chemical species, often referred to as a substrate. The very simplest version (which 




(1.26) 




(1.27) 
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is not actually cooperative!) is the Michaelis-Menten reaction scheme 



S + E^C%E + P, (1.28) 

k—i 

which represents the reaction of a substrate (S) with an enzyme (E) to form a complex 
(C) which produces a product (P), together with the original enzyme. The enzyme 
thus catalyses the reaction (i.e., it is not absorbed by the reaction). 

Although the law of mass action is applied to these reactions separately, it is 
usually the case that the first reaction (the formation of the complex) is very fast, 
essentially because there is normally not much enzyme present. As a consequence 
of this, the forward rate k\SE is approximately equal to the backward rate k_\C. 
In addition, the total enzyme is conserved, either as free enzyme E or in the bound 
form C, so that E + C = E is constant. It follows from these observations that the 
product rate formation, i. e., the rate of reaction, is 

where K\ = k-\/k\. This provides a reaction- based explanation for saturating rate 
kinetics. 

The actual phenomenon of cooperativity refers to the ability of some enzymes to 
have multiple binding sites for certain substrates. For example, with two binding 
sites we might consider the reaction scheme 

S + E^d-tE + P, 
k-i 

S + d 5= C 2 % C 2 + P. (1.30) 

k-3 

It is easy to analyse this in the same way, assuming the binding reactions are fast. 
We find that the overall production rate is 

_ (^K.S + hK.KsS^Eo 

where K\ = fc-i/fci and K 3 = k_ 3 /k 3 . The predation rate in (1.26) is a particular 
case of this, if K 3 is large enough, 

It is easy to extend these ideas to enzymes with n binding sites, and we then 
find that the reaction rates are expressible as ratios of n-th degree polynomials. In 
particular, such models can produce the rates described by Hill functions, in the form 

r * FT*' (133) 

such models are often used as empirical rate terms, for example in respiratory control 
modelling. 
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1.5 Notes and references 



There are quite a lot of books about mathematical modelling, but they do vary in 
range. What you certainly don't want are books on simple methods: how to do the 
Laplace transform, and so on. There are a number of relatively recent books on the 
principles of mathematical modelling, for example Tayler (1986), Haberman (1998), 
Fowkes and Mahony (1994), Howison (2005); these are aimed at undergraduates, as 
is the present text. Fowler (1997) aims at a more graduate level. The classic of the 
type is the book by Lin and Segel (1974), now very dated but certainly an absolute 
model of clarity. 

There are obviously many books which describe the technical material which is 
used here, but not so many which are written from the viewpoint of the practising 
applied mathematician. Of several fitting this latter description, we mention Keener 
(2000), Ockendon et al. 2003, and Carrier and Pearson (1976). 



The SIR model was introduced to epidemiology by Kermack and McKendrick (1927). 
It is described in the book by Murray (2002), where various extensions to the model 
are described. Recovery of the susceptible population (for example by population 
growth) can lead to oscillations, indicative of repetitive outbreaks. If these occur in 
a spatially distributed domain, then population migration (modelled as a diffusion 
term) can cause an epidemic to propagate as a travelling wave. Apparently this is 
how the Black Death spread in Europe in the fourteenth century. 

The short book by Hoppensteadt (1975) describes age dependent population mod- 
els, as well as the SIR model. 

1.5.2 Spruce bud worm model 

Murray's (2002) book also describes the simplest version of the spruce budworm 
model. This itself is a simplification of a three variable continuous model described 
by Ludwig et al. (1979), which itself is an applied mathematician's attempt to propose 
a simpler version of the complex simulation model of Jones (1979). The Ludwig model 
is analysed in greater depth by Fowler (1997), and the Jones 'budworm site model' 
has been analysed by Hassell et al. (1999). Royama (1992) discusses the ecological 
background of the problem. 

Exercises 

1.1 Consider the discrete population model (1.1): 



1.5.1 SIR model 




— hu. 
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By writing u n = Lw n for some suitable choice of L, show that the model takes 
the form 

w n+ i = Xw n (l - w n ), 

and determine A. What is the effect of increasing h on the behaviour of the 
population? What happens if h > r — 1? 

1.2 Consider the continuous population model (1.2): 

du I ' u 



= pu (l 



dt "~V CJ ^ 

By writing u = Mw for some suitable choice of M, show that the model takes 
the form 

dw 

- = aw(l-w), 

and determine a. What is the effective of increasing n on the population? What 
happens if p < /i? 

1.3 Starting from an integral conservation law, derive the heat equation in the form 

^(pc p T) + V.[pc p Tu] = V.[kVT] + Q, 
and the mass conservation equation in the form 

! + V.(pu) = 0. 

1.4 The SIR model of an epidemic is 

S = -kSI, 
I = kSI-rl, 
R = rl, 

where S represents the susceptible population, / the infected population, and 
R the recovered (or removed) population. Suppose that at t — 0, S — So, I is 
very small (perhaps a single individual), and R = 0. Show that an epidemic 
will occur if the reproductive rate 

R= k_So 
r 

is greater than one, and find / as a function of S. 

1.5 A model for the reaction of a substrate with an enzyme with two binding sites 
is given by 

fc i k 2 

S + E^d-tE + P, 

k-i 
k 3 k 4 

S + d ^ C 2 C 2 + P 

k- 3 
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If the reversible reaction rates k±i and k± 3 are very fast, show that the rate of 
formation of product P is approximately 

1 + K.S + K^S 2 ' 
where Ki = k_\jk\ and K 3 = k_ 3 /k 3 . 
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Chapter 2 



Non-dimensionalisation and 
approximation 



Once we have a model, we have to try and solve it. There are two kinds of solutions: 
exact, analytical solutions, and approximate solutions. Exact solutions are explicit 
formulae; for example we can exactly solve quadratic equations, and certain differen- 
tial equations, such as that describing simple harmonic motion. We also consider that 
solutions such as Taylor series constitute analytic solutions: they can be computed 
to arbitrary accuracy. The same applies to quadratures, such as the solution of 



Approximate solutions are those where one solves an approximate equation, or an 
approximating sequence of equations. Approximate methods are best applied when 
the approximation is based on the size of certain terms. In this chapter we will 
illustrate the use of such methods, firstly on simple algebraic equations, and then on 
some differential equations. The whole subject of perturbation theory is extensive, 
and a thorough discussion is beyond the scope of these notes. 

2.1 Non-dimensionalisation 

In order to approximate a solution, we need to be able to neglect terms which are 
small. This raises a concept of fundamental importance, which is that 'small' and 
'large' are adjectives which can only apply quantitatively when a comparison is made 
between quantities of the same dimension. An equivalent assertion is that we can 
only make approximations based on the small size of parameters if these parameters 
are dimensionless. It makes no intrinsic sense to say a quantity is small if it still has 
dimensions. A speed of 1 cm s _1 is small if you are a rocket, but large if you are a 



du 
lit 



f(u), u(0)=u , 



(2.1) 



which has an implicitly defined solution 




(2.2) 
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giant African land snail. An ice sheet is thick if you are a human being, but thin if 
you are a planet. So we always associate large or small quantities with dimensionless 
numbers, that is, parameters which measure quantities with respect to some reference 
value. The process of isolating these parameters is called non-dimensionalisation. 



2.1.1 The wave equation 

Putting a mathematical model into non-dimensional form is fundamental. Although 
technically trivial, there is a certain art to the process of non-dimensionalisation, and 
the associated concept of scaling, and the only real way to learn how to do it is by 
example. Let us begin with a simple model, the wave equation, which one learns how 
to derive in first year applied mathematics courses. We suppose a string, for example 
a guitar string, is wound between two points a distance I apart. If the tension in the 
string is T and its density (per unit length) is p, then an application of Newton's 
second law to an infinitesimal segment of the string leads to the equation 

d 2 y d 2 y 

p W =T d^> y = on X = (U (2 - 3) 

where x is distance along the string, and y is its transverse displacement. 

The main assumption that is usually stated in deriving this equation is that the 
displacement y is small. However, there are at least two other implicit assumptions 
which are made. One is that gravity is unimportant; the other is that there is no 
longitudinal displacement. 

For a guitar string, these seem to be reasonable assumptions, but why? We expect 
the effect of gravity to be a deviation of the displacement from the vertical, and this 
is evidently valid for the guitar string. It is not valid for the hanging chain, or for 
the wire between telegraph poles. Why? I would say, for the chain, the density is 
too large; for the telegraph wire, the distance I is too large; while the guitar string 
is straight because it is tight: T is large. These facts suggest that the 'size' of the 
gravitational acceleration g may in fact be measured by the dimensionless parameter 
pgl/T, which appears to be the only independent dimensionless parameter which can 
be formed from those in the model if we include gravity. 

How can this suspicion be confirmed? From first principles, we derive the wave 
equation, including gravity, in the form 

d 2 y d 2 y 

p W = T dl 2 ~ p9 ' y = on X = 0J - (2 - 4) 

Next we write the model in dimensionless form. We do this by non-dimensionalising 
the variables so that they are numerically of order one (written 0(1)). Specifically, 
we write 

x = lx\ y = y y\ t = (l/c)t\ (2.5) 



where 

It 



(2.6) 
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2-^*2 ft (2-7) 



is the wave speed, and y is a measure of the displacement: for example, it could be 
the maximum initial displacement. The dimensionless model is then obtained in the 
form 

d 2 y* d 2 y 

dt* 2 ~ dx 
where 

Ty 

with y* = 0(1) initially, and y* = on x = 0, 1. All of the terms in the equations 
and in the initial and boundary conditions are dimensionless, and all the coefficients 
which appear (such as (5) are dimensionless. 

It is conventional at this dimensionless stage to dispense with the asterisks in 
writing the variables, and this we now do. The process of choosing particular scales 
for the variables, or scaling, is motivated by the following ideas. Firstly, there is a 
natural length scale to the problem, /, which is the dimension of the geometric domain 
on which the problem is to be solved. Further, there is a natural length scale for the 
displacement, yo, which is present in the initial conditions. Next, of the three terms 
in (2.4), we anticipate (at least for the guitar string) that the gravity term will be 
'small'. It follows that the other two should be the same size, and we choose the 
time scale so that these two terms 'balance', that is their dimensionless scales (here 
Ty 2 /l 2 ) are the same. When the model is written in the dimensionless form (2.7), we 
then have equal dimensionless coefficients multiplying these two terms: here they are 
both equal to one. 

The final, essential idea is that, in general, a dimensionless function u(x, t) which 
varies by 0(1) over an x range of 0(1) will have derivatives of 0(1). This is true, 
for example, for sinx, e~ x , and x 2 : it is not true for the function e~ 10x , which varies 
rapidly over a distance of order x ~ 0.1 near x = 0. With this assumption, the 
derivative terms d 2 y/dt 2 and d 2 y/dx 2 are 0(1), and it follows that the relative size 
of the gravity term is given by (5. Thus gravity is negligible if j5 <C 1, and indeed this 
means large tension, small density or short length, as we surmised. 

2.1.2 The heat equation 

Next we consider a form of the heat equation, (1.8). We write it in the form (assuming 
density p and specific heat c p are constant) 

dT 

— + u.VT = kV 2 T + H, (2.9) 
dt 

where H = Qj pc p . We have assumed V.u = 0, which follows from the conservation 
of mass in the form 

^ + V.(pu) = 0, (2.10) 

together with p = constant. 

Suppose we are to solve (2.9) in a domain D of linear dimension /, on the boundary 
of which we prescribe 

T = T B on dD, (2.11) 
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where Tg is constant. We also have an initial condition 

T = T (x) in D, t = 0, (2.12) 

and we suppose u is given, of order U. 

We can make the variables dimensionless in the following way: 

x = Zx*, t = [t]t', T = T B + (AT)T*. (2.13) 

Again, we do this in order that both dependent and independent variables be of 
numerical order one, 0(1). If we can do this, then we would suppose a priori that 
derivatives such as V*T* (V = / _1 V*) will also be of numerical 0(1), and the size 
of various terms will be reflected in the dimensionless parameters which occur. 

In writing (2.13), it is clear that I is a suitable length scale, as it is the size of D. 
For example, if D was a sphere we might take I as its radius or diameter. We also 
suppose that the origin is in D; if not, we could write x = xo + Zx*, where xo € D: 
evidently x* = 0(1) in D. 

A similar motivation underlies the choice of an 'origin shift' for T. In the absence 
of a heat source, the temperature will tend to the uniform state T = Tg as t — > oo. 
If H ^ 0, the final state will be raised above Tg (if H > 0) by an amount dependent 
on H. We take AT to represent this amount, but we do not know what it is in 
advance — we will choose it by scaling. The subtraction of Tb from T before non- 
dimensionalisation is because the model for T contains only derivatives of T, so that 
it is really the variation of T about Tb which we wish to scale. 

In a similar way, the time scale [t] is not prescribed in advance, and we will choose 
it also by scaling, in due course. 

With the substitutions in (2.13), the heat equation (2.9) can be written in the 
form 

/ I 2 \ dT* (Ul\ , „, m , „,„ , / HI 2 \ . . 

U j W + \ k J = v* 2 r + , (2.i4) 

where we have written u = Uu*, so that u* = 0(1). This equation is dimensionless, 
and the bracketed parameters are dimensionless. They are somewhat arbitrary, since 
[t] and AT have not yet been chosen: we now do so by scaling. 

The solution of the equation can depend only on the dimensionless parameters. It 
is thus convenient to choose [t] and AT so that two of these are set to some convenient 
value. There is no unique way to do this. 

The temperature scale AT appears only in the source term. Since it is this which 
determines the temperature rise, it is natural to choose 

HI 2 

AT = — . (2.15) 

K 

It is also customary to choose the time scale so that the two terms of the advective 
derivative on the left of (2.14) are the same size, and this gives the convective time 
scale 

[t] = 1 (2.16) 
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Finally, we remove the asterisks. When this is done, the dimensionless equation takes 
the form 



Pe 

where the Peclet number is 



w +uVT 



= V 2 T + 1, (2.17) 



Pe = — , (2.18) 

K 

and the solution of the model depends only on this parameter (as well as the initial 
condition). The boundary condition is 

T = on 3D, (2.19) 



T = 0(x) at t = 0, (2.20) 



,(x) = T «^. (2.2!) 



and the initial condition is 
where 

2.2 Scaling 

A well-scaled problem generally refers to a model in which the dimensionless param- 
eters are 0(1) or less. Evidently, this can be ensured simply by dividing through 
by the largest parameter in any equation. More importantly, if parameters are nu- 
merically small, then (as we discuss below) approximate solutions can be obtained 
by neglecting them. The problem is well-scaled if the resulting approximation makes 
sense. For example, (2.17) is well-scaled for any value of Pe. However, the problem 
eT t = eV 2 T + 1, with e <C 1, is not well scaled. One makes a problem well-scaled in 
this situation by rescaling the variables, and we now consider the wave equation (2.7) 
again in this light. 

2.2.1 The wave equation, again 

In the statement of the wave equation with gravity, there are in fact two dimensionless 
parameters: 

B = f, e=*. (2.22) 

The parameter e is a measure of the amplitude of the motion, and it is on the basis 
that £<1 that we derive the linear wave equation in the first place. The parameter 
/3 — B/e, so the assumption of negligble gravity is equaivalent to the assumption that 

5<£<1. 

If B ~ e, then /3 = 0(1). The problem is sensibly scaled, but gravity is no longer 
negligible. There is a steady state y — \fix(l — x) (the hanging chain), and, because 
the equation is linear, the string simply oscillates about this steady state. 

Now suppose that (3 ^> 1. The model is now not correctly scaled (because the 
limit (3 — > oo gives no sensible approximation). In fact the model suggests that a 
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steady state will have y ~ j5, and that the string will oscillate about this steady state 
with a similar amplitude. In order to obtain a sensibly scaled problem, we simply 
have to rescale the displacement by writing y = (3Y, and the discussion can proceed 
as before, except that the initial condition has Y<1. 

It seems that all is now well; we have discussed the cases f3 <C 1, f3 — 0(1), 
and |9 > 1. However, ther is more concern when f3 becomes as large as 0(l/e). 
In this case, the model suggest oscillations about a steady state of order j5 ~ 1/e, 
and dimensionally of order /. Because of this, the basis of the original derivation is 
suspect, and the model must be re-derived: in fact this can be done (see question 
2.1), so that the model equation (2.7) remains valid. However, the initial value scale 
yo is not now an appropriate scale for y\ the appropriate (dimensional) scale is y ~ 
and this a posteriori adjustment is the rescaling alluded to above. In dimensionless 
terms, we rescale the model by writing y = (l/y )Y = Y/e, and we find 

d 2 Y d 2 Y 

and this version of the model is appropriate if B = 0(1). If B 3> 1, a further rescaling 
simply takes Y ~ B. 



2.3 Simple approximation methods 

Suppose we wish to solve the equation 

x 3 - ex -1 = 0, (2.24) 

where e is relatively small: for example if £ = 0.1. Formulae do in fact exist for 
writing solutions of cubic (and also quartic) equations, but they are fairly unpleasant 
and are rarely used. A much better way is to use an approximation method, based 
on the idea that the parameter e in equation (2.24) is small. 

Graphically, it is clear (see figure 2.1) that when e is small, (2.24) will have a 
unique, positive real root, and in fact it will be close to x = 1 (since 1 + ex « 1). 
(In passing, note that there will be exactly one root for e < e c , where e c is the value 
corresponding to tangency of the line with the curve; at this value x 3 = 1 + ex and 
3x 2 = £, so that e c = 3/2 2/3 w 1.89.) 

A simple iterative method to solve (2.24) is 

x n+1 = {l + ex n ) 113 . (2.25) 

If we choose e = 0.1 and xo = 1, then successively 

xi = 1.0322801, 

x 2 = 1.0332889, 

x 3 = 1.0333204, 

x A = 1.0333214, (2.26) 
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Figure 2.1: The graphs of x 3 and 1 + ex, with e = 0.1. 

and this last value is the root. It is of course trivial to compute the positive root for a 
range of e, but it would be convenient to have an analytic (as opposed to numerical) 
approximation. We construct this by using a perturbation method. 
We can use the binomial expansion to write (2.24) in the form 

x = (1 + ex) 1 ' 3 = 1 + \ex - \e 2 x 2 + ^eV .... (2.27) 

Since e <C 1, we see that isil [+0(e): that is, terms of order e, i.e. of size about e]. 
A better estimation is then 

x « 1 + \ex « 1 + \e, (2.28) 

and we can see that this relatively crude approximation is in fact accurate to four 
decimal places when e = 0.1! Repeating this idea, we would have 

/y> r^J 1 _|_ 1 <-> rp 1 r-2 /T ,2 

Jb r-^j ± -\- -fcX — — c JU 

W 1 + |£(1+|£)-|£ 2 (1 + ...) 2 ^1 + |£ (2.29) 

(there is no 0(e 2 ) term), but a more methodical procedure is to anticipate (by in- 
spection) that the root can be written in the form of a series 

x = x + ex\ + e 2 x 2 + . . . ; (2.30) 

we substitute this into (2.24), expand in powers of e, and equate coefficients of powers 
of e: a little thought indicates why this procedure is necessary. For (2.24), we thus 
have 

(x + exi + e 2 x 2 + . . .) 3 - + ex l + ...)- 1 = 0, (2.31) 
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whence 



(xl — 1) + e(3xlxi — xo) + e 2 (3xlx 2 + 3x x 2 — xi) 

+ e 3 (3xlx 3 + 3x x 1 x 2 +x 3 1 -x 2 ) + ... = 0, (2.32) 



so that, sequentially, 



4-i = o, 

3xfai — x = 0, 
3x^x 2 + 3xox 2 — xi = 0, 
3xgX 3 + 3x xix 2 + Xi - x 2 = 0, 



from which we obtain 



and hence the root is 



x 



(2.33) 



x = 1, xi = |, x 2 =0, x 3 = (2.34) 



1 + |e- ^£ 3 + 0(£ 4 ). (2.35) 

With £ = 0.1, we have x pa 1.033321: practically, the exact result. Even, for e — 1, 
the approximate result is 1.321, while the exact root is 1.3247. 



Singular approximations 

The approximation above is called a regular approximation, because the limit when 
£ = gives an approximation to the root. Now consider the cubic 

£x 3 - x - 1 = 0, e < 1. (2.36) 

Graphically (figure 2.2), there are clearly three real roots. One of these is near —1 
and can be recovered by a regular approximation. The others are at large values of 
|x|, and are determined by balancing ex 3 with x; that is, £x 3 ~ x when x ~ £ -1 / 2 , so 
we first write 

x = e-^X, (2.37) 



and then 

X 3 - X - £ 1/2 = 0, (2.38) 

with approximate roots X pa 0, ±1. The root X pa corresponds to x pa —1 and is 
determined by the regular approximation for x. The larger roots are determined by 
a regular approximation of (2.38), as a power series in £ 1//2 . Thus 

X = Xo + e 1/2 X 1 +eX 2 + ..., (2.39) 

and substituting this into (2.38) and equating powers of £ 1//2 , leads to the approximate 
solutions (written in terms of x) 

x ~ ± \ + |t|v / £+ 3 ^£ + ---- (2.40) 
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Figure 2.2: The graph of ex 3 — x — 1, for e — 0.1. 



The upper sign gives the positive root, the lower the negative root. For example if 
e = 0.1, the approximate positive root from (2.40) is 3.565567, while the exact root is 
3.5770894. For e = 1, the approximate root is 1.1468753, while as we have seen, the 
exact root is 1.3247. This is less good, but can be improved by taking further terms 
in the series (2.40). We see that approximation methods can provide a very useful 
way of solving algebraic equations. 
Now suppose we wish to solve 

tanx = tanhx. (2-41) 

Each function is odd, so x = is a solution, and graphically (figure 2.3) it is clear 
that there is a sequence of positive roots x\,x\,... (and thus also negative roots 
— x\, —x\, ■■■), and that as n — > 00, 

<«(ra + J)7r. (2.42) 



Suppose we put 
then 



so (2.41) is 



x' n = (n+\)n + 6; (2.43) 



(-1)" 

sinx* = — y= (cosfl + sinfl), 

cos< = ^^(cosfl-sinfl), (2.44) 
V2 



cos + sin 6 1-e 2x * n . „ . 

tan< = — — — _ 2x , . 2.45 

cos — smu 1 + e 2x n 
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Figure 2.3: The graphs of tan a: and tanhx. 



We expect 6<1, and also x* n 3> 1, and thus, since 

cos# ~ i-|# 2 , 

sin0 ~ 6, (2.46) 
we have, with use of the binomial expansion for (1 + 

]:^;;;^;;; = a - e - 2< )a - e - 2 < + . . .) 

(1 + - ^ 2 . . .)(1 + {0 + \6 2 + . . .} + {<9 2 + ...}...) = 1 - 2e" 2 < . . . 

=> 1 + 29 + 29 2 + . . . = 1 - 2e- {2n+1 ^- 2d . . . 

^ + ^ 2 + ... = _ e -( 2 "+D 7r (l-20...), (2.47) 

so that, finally, 

~ - e "( 2n+ ^, (2.48) 

and 

< ~ ( n + J)t - e"( 2n+ 5)\ (2.49) 
Numerical approximation 

Iterative numerical methods are widely used to solve algebraic equations. A general 
iterative method to find a root of f(x) = is to define a sequence x — xq, xi, . . . , 
satisfying x n+ i — g(x n ), where the function g is chosen so that x = g(x) if f(x) = 0: 
for example, g{x) = fix) + x. A simple iterative method to solve 

L(x) = R(x) (2.50) 

is as follows: define 

L(x r+1 ) = R(x r ), (2.51) 
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i.e. x r+ i = Lr 1 o R(x r ). The sequence will converge if \{L _1 o R)'\ < 1 at a root, and 
using the chain rule (via / = L _1 o R if L[f(x)] = R(x), so L'f = R 1 , and at a root 
f(x) — x, f — R'/L 1 ) we find that this is \L'\ > \R'\. Consulting figure 2.3, we see 
that the iteration 

x r+ i = tan _1 [tanha: r ] + nir (2.52) 

will converge to x* n (since tan -1 is defined to be less than ir/2 on a calculator). 
The lowest root (n = 1) is approximated by 

x\ « |tt - e~ 57r/2 = 3.9266026 (2.53) 

compared with the exact value 3.9266024. Here an approximation based on the limit 
n — > oo is accurate even when n = 1. 



2.4 Perturbation methods 

Let us consider (2.17) with (2.19) and (2.20), and suppose that O < 0(1). If Pe -C 1, 
we obtain an approximation by putting Pe = 0: V 2 T + 1 w 0. Evidently, we cannot 
satisfy the initial condition, and this suggests that we rescale t: put t = Per, so that 
(approximately) 

dT 

— = V 2 T + 1; (2.54) 

OT 

now we can satisfy the initial condition (at r = 0) too. Often one abbreviates the 
rescaling by simply saying, 'rescale t ~ Pe, so that T t w V 2 T + 1'. 




Figure 2.4: Sub-characteristics and boundary layer for the equation (2.17). The sub- 
characteristics are the flow lines cht/dt = u, and the boundary layer (of thickness 
0(1/ Pe)) is on the part of the boundary where the flow lines terminate. 

On the other hand, if Pe ^> 1, then T t + u.VT w 0, and we can satisfy the initial 
condition but not the boundary condition on all of dD, since the approximating 
equation is hyperbolic (its characteristics are called 'sub-characteristics'). To remedy 
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this, one has to rescale x near the part of the boundary where the boundary condition 
is not satisfied. This gives a spatially thin region, called (evidently) a boundary layer, 
of thickness 1/Pe (see figure 2.4). 

The other possibility is if ^> 1, say 9 ~ 6q S> 1. We discuss only the case Pe ^> 1 
(see also exercise 1.2). Since T ~ 6q initially, we need to rescale T, say T = OqT. 
then Pe[f t + u.VT] = V 2 T + Bq 1 , and with f = 0(1), we have f t + u.VT w 
for Pe 3> 1. The initial function is simply advected along the flow lines (sub- 
characteristics), and the boundary condition T = is advected across D. In a 
time of 0(1), the initial condition is 'washed out' of the domain. Following this, 
we revert to T, thus T t + u.VT = Pe _1 (V 2 T + 1). Evidently T will remain w 
in most of D, with a boundary layer near the boundary as shown in figure 2.4. If 
n is the coordinate normal to 3D in this layer, then u.VT ~ u n dT/dn ~ PeT, 
Pe _1 V 2 T ~ Pe~ 1 d 2 T/dn 2 ~ PeT, and in the steady state, these must balance the 
source term Pe -1 : thus in fact, the final state has the rescaled T ~ Pe -2 , and this 
applies also for 6 < 0(1). 

These ideas of perturbation methods are very powerful, but a full exposition is 
beyond the scope of these notes. Nevertheless, they will relentlessly inform our dis- 
cussion. While it is possible to use formal perturbation expansions, it is sufficient in 
many cases to give more heuristic forms of argument, and this will typically be the 
style we choose. 

2.5 Notes and references 

Non-dimensionalisation is technically trivial, analytically essential, and difficult to get 
right. You won't find many books that dwell on it, because of its apparent triviality. 
But, it is also clear that you can get it very wrong. The clearest precise exposition 
for the uninitiated is in the book by Lin and Segel (1974), and other than that it is 
largely a matter of example and practice. What is clear is that you can easily get it 
wrong! 



Exercises 

2.1 Derive the wave equation describing oscillations of a string of length I from first 
principles, when gravity is included, assuming displacements are small and of 
order y . Show how to non-dimensionalise the equation to obtain the form 

d 2 y = d 2 y 
dt 2 dx 2 P ' 

and define (3. 

Now suppose that yo ~ I. Suppose that in the unstretched state where the 
displacement y = 0, the density po is constant. By careful consideration of 
the application of Newton's second law to an infinitesimal segment of length 
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ds (stretched from its original length dx), show that the (dimensional) wave 
equation can be derived in the form 



d 2 y m d 2 y 



dt 2 



T °dx~2-P°9, 



assuming only that displacements are purely vertical, where To is the horizontal 
component of the tension exerted at the fixed end points at x = 0, x = I. 

Non-dimensionalise the model in this case and describe the form of the resulting 
oscillation. 

2.2 A population of size N is subject to immigration at rate /, and mutual pair 
destruction at a rate kN 2 , so that N — I — kN 2 . By appropriate scaling of the 



variables, show that the model can be written in the form x 



x . 



2.3 A model describing the interactions of populations of a herbivore (H) and plank- 
ton (P) is 



dP 

~dt 
dH 

~dt 



= rP 
= DH 



(K-P)- 
P 



BE 



C + P 



C + P 
AH 



Explain what the terms in the equations mean, and suggest ecological reasons 
for their form. 

Show how to non-dimensionalise the model to obtain the dimensionless pair of 
equations 



x = x 



y 



ay 



k — x 

x 



y 



Li + x 



l + x 
- ay 



and define the dimensionless parameters k, a amd a. By identifying the units 
of the physical constants in the model, show that these parameters are indeed 
dimensionless. 



2.4 Suppose 



with 



Pe 



dT 
~dt 



+ u.VT 



V 2 T + 1 in D, 



T = on 3D, 

T = # o 0(x) in D at t = 0, 

and © = 0(1), #o 3> 1, -Pe -C 1. Discuss appropriate scales for the various 
phases of the solution. 
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2.5 Explain why the iterative method 

x n+1 = (1 + ex n ) 1/3 

used to solve (2.24) will converge to its solution. Does this depend on the value 
ofe? 

2.6 (i) Find approximations to the solution of 

ex 3 - x - 1 = 0, £«1, 

which is close to x = — 1. Compare with the numerical solution when e = 0.1; 
e = 0.01. 

(ii) Use perturbation methods to find approximate roots to the equation 

xe~ x = e, < £ « 1. 

(Use graphical methods to find the location of the roots. For the larger root, 
take logs and note that if x ^> 1, then x ^> lnx.) 



2.7 Each of the equations 



Z 5 - £Z - 1 = 0, 

ez 5 - z - 1 = 0, 



has five (possibly complex) roots. Find approximations to these if e <C 1. Can 
you refine the approximations? 

2.8 Suppose that m satisfies the polynomial 

— m 3 + e(A + B)m 2 - (A - A)m - 2 7 (A + B) = 0, (*) 

7 

where a, 7, A, B, A are 0(1) and £<Cl. Show that the roots are approximately 
given by 



m ~ , m± w ±i ( — — 

A — A \ sol 



1/2 



Now suppose that 7 = g/e and m = Mje. Show that M satisfies 

-M 3 + (A + B)M 2 + {A- X)M - 2g(\ + B) = 0. (**) 

Suppose that g is large. Show that the roots are approximately given in this 
case by 

(A + B)g 



M± fa ±J2g, M « 

v a 

By consideration of the roots of (**) as g decreases, can you ascertain the path 
of the roots in the complex plane, and how M± and M are transformed to the 
roots m± and mo of (*)? 
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2.9 The Lorenz equations axe 

x = —ax + ay, 
y = (r - z)x - y, 
z = xy — bz. 

These equations are already dimensionless, and r, b and a are dimensionless 
parameters. Is it possible to rescale the equations in order to reduce the number 
of independent parameters? If so, how? If not, why not? 

Suppose the parameter r is large. Show how to rescale the equations so that 
the damping terms in the equations are small (of size e = 1/ ^/r), and show that 
if the terms in e are put to zero, there are two first integrals of the equations. 
What does this imply about the solutions when e = 0? 

2.10 The SIR model of an epidemic of a fatal disease is 

S = -kSI, 
I = kSI-rl, 
R = rl, 

where S represents the susceptible population, / the infected population, and 
R the removed (dead) population. Suppose that at t = 0, S = So, I is very 
small (perhaps a single individual), and R = 0. 

Non-dimensionalise the model to obtain the form 

s = sf, 



f 

where the dimensionless parameter 



R= kSo 
r 

is called the reproductive rate. 

Find a first integral for / in terms of s, and show that an epidemic will occur if 
R > 1. 

Now suppose that the epidemic is weak, in the sense that R — 1 = e <C 1. Show 
that a suitable approximation of the model can be made by defining 

T 

s = 1 - ea, f = e 2 F, t = — , 

£ 

and show that the (scaled) death rate F of the resultant approximation is ap- 
proximately given by 

F = |sech 2 |T. 
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Chapter 3 

Perturbation methods 



In chapter 2 we were introduced to the basic ideas of approximation techniques, 
which are central to the art of the practical applied mathematician. In this chapter 
we introduce some of the formal analytic methods which have been devised to provide 
approximate solutions to problems which we will suppose are framed as differential 
equations. We begin with some generalities, but the subject is best taught and learned 
by way of example, because by its nature there is no general theory which will provides 
a method of solution for all problems. 

Suppose we have an ordinary differential equation for a function u(t), which we 
write in the formal way 

N(D,e;u) =0; (3.1) 
here N denotes a nonlinear differential operator of the ordinary derivative D = 

at 

Often N will be a polynomial function of both D and e. If we suppose that N is 
analytic in e, then we can expand the equation as 

N(D,e;u) = N {D ;u) + eN x (D ;u) + . . . = 0, (3.2) 

and this suggests that we seek a perturbation series for the solution of the form 



u 



X>X, (3.3) 



s=0 



and if we substitute this into (3.2), then we find that 

Y,e r+a N ri8 (D;u a ,u a - 1 ,...)=0. (3.4) 

r,s 

Here we have expanded each operator term as 

(oo \ oo 

= Y, e 8 N r , a {D ;u s ,u s . u . . .). (3.5) 

Note that each operator N r ^ s depends on functions with indices up to s. By 
changing the indices of summation, we can write (3.4) in the form 

oo 

J2e j N j (D;u j ,u j _ 1 ,...)=0, (3.6) 

3=0 
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where 

j 

Nj(D;Uj,Uj_ u ...) = Nj- SjS (D ;u a ,u a -i, . . .)■ (3.7) 

The formal procedure is now to equate powers of e to zero, yielding a sequence of 
equations 

N j (D;u j ,u j - 1 ,...) = 0. (3.8) 
It is clear that we can, at least in principle, solve these sequentially for Uj to find 

U j = ^(«i-l,«i-2, • • •)• (3-9) 

This sequential procedure then yields our solution. 

Apart from the fact that there are differential operators involved, what we have 
described above is simply an iterative procedure for evaluating a Taylor series solution 
for u in terms of e. What makes perturbation methods interesting is that the series 
that is so constructed is rarely convergent. The more usual divergent power series is 
called an asymptotic series, or an asymptotic expansion, and is written 

oo 

u~J2 £Su s- ( 3 - 10 ) 

The symbol '~' denotes asymptotic equivalence, and in general we write u ~ v (as 
£ — > 0, for example) if the ratio u/v tends to a fixed limit as e — > 0; we also write 
u = 0(v) in this case. Other notation is to write u = o(v) (or u <C v ) if u/v — > as 
£ -)• 0. 

Despite its divergence, the series (3.10) can still yield a good approximation to 
the solution when e is small. The reason for this is that the terms will then initially 

N 

decrease, and so the finite sum ^ u s will yield an increasingly accurate solution if iV 

o 

is fixed and e is decreased. 

Why should the series for u diverge if the original operator N is analytic in u and 
£? The divergence of an asymptotic series such as (3.10) indicates that the function 
u(t, e) has a singularity at e = 0. This must certainly be the case if u does not 
approach the solution uq of the unperturbed equation 

N (D;u ) = 0, (3.11) 

and this is often the case. It is, for example, generally true for singularly perturbed 
problems, and these are precisely the types of problem which are of most interest 
(because they are the most difficult). 

3.1 Regular perturbation theory 

We illustrate some of the above abstract discussion by consideration of some specific 
examples. Let us consider the boundary value problem 

u" - u = eu 2 , «(0) = 0, u{l) = 1, (3.12) 
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where £<1, and we will denote the independent variable as x, since boundary value 
problems usually arise in a spatial context. 

Following our description above, we seek a solution in the form (3.10), and sub- 
stituting this into (3.12) and equating powers of e, we find the sequence of equations 

u'q -u = 0, 

U'l -Ul = Uq, 

u" 2 — u 2 = 2u ui, (3.13) 

and so on. The boundary conditions are also expanded in this way, so that we have 

uo = u i = u 2 = at x = 0, 
uo = 1, «i = «2 = at x — 1. (3-14) 

It is easy, if messy, to solve this sequence of equations. We have 
sinhx 

Uq 



sinh 1 



-\ + \ cosh 2x + \ cosh x (| - | cosh2 - | cosh l) sinh x 

1,1 = rfi + ss^i • (ai5) 

for example. We can see that «o contains terms in e x and e~ x , while u\ contains 
terms from e 2x down to e~ 2x . It is not difficult to show by induction that u r contains 
terms from e^ r+1 ^ down to e~( r+1 K This suggests that the perturbation series will 
be convergent for sufficiently small e; see also question 3.2. Note that this will not 
be true for an infinite interval of integration, for example if we prescribe two initial 
conditions at x = and require solutions valid in < x < oo. We discuss this case 
further in section 3.3 below. 



3.2 Singular perturbation theory: boundary layers 

Next we consider an equation which is fairly similar to (3.12), but with the small 
coefficient e multiplying the highest derivative: 

eu" - u' = u 2 , u(0) = 0, u(l) = 1, (3.16) 

and £< 1. 

We use the same procedure as before, writing 

u ~ u (x) + eu 1 (x) + . . . , (3.17) 

and then equating terms of the same order to obtain 

u'q + uI = 0, 
u[ + 2m mi = Uq, 

u' 2 + 2uqU2 = u'[ — ul, (3.18) 
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and so on. The associated boundary conditions are 



u = Ui = u 2 = at x = 0, 
uo = 1, «i = «2 = at x = 1. (3.19) 

The essential distinguishing feature about this problem is that the small coefficient 
multiplies the highest derivative of the equation. Because of this, the leading order 
problem (for m ) is of lower degree than the full equation, and in general this means 
that it is not possible to satisfy all of the boundary conditions. 

An issue of choice then arises, because this leading order 'outer' solution may be 
chosen differently depending on which boundary conditions are satisfied. In addition, 
nonlinearity of the equation may allow the possibility of genuinely multiple solutions. 

In the present case, we have at leading order 

u' = -u 2 . (3.20) 

Let us suppose, apparently arbitrarily, that we apply the left hand boundary condition 

u = at x = 0, (3.21) 

and in general we will then aim to satisfy the boundary conditions on u n at x = 0. 
As we shall see, the issue of choice is generally resolved by requiring consistency of 
the solution procedure. 

Picard's theorem tells us that 

u = (3.22) 

is the only solution of (3.20) with (3.21). Since uo is constant, it is not difficult to 
see that the successive terms u\, u 2 and the rest will also all be zero, so that this 
outer solution is given by u = to all orders of e. Evidently, u = cannot be 
the exact solution (though it exactly solves the equation) because it does not satisfy 
the boundary condition at x = 1. Here we have an early glimpse of the notion of 
exponential asymptotics: the outer solution for u is smaller than any algebraic power 

of e, and this suggests that it may be exponentially small, or u = exp — O ( - 

L Ve, _ 

In any event, we have an approximation to the solution which does not satisfy the 
boundary condition at x = 1. This is associated with the loss of the highest derivative, 
and such problems are called singular perturbation problems, presumably because 
the resulting perturbation solution is singular (having, in this discontinuity at 

x = 1). In order to resolve the failure of the method, we need to be able to bring back 
the highest derivative term into the approximation scheme, and the way in which we 
can formally do this is by adjusting the length scale. We do this near x — 1, since 
that this where we have the problem. We put 

x = l-6X, (3.23) 

and suppose that 5, as yet unidentified, is small. Changing the independent variable 
in (3.16), we find that u, now taken as a function of X, satisfies the equation 

j 2 u" + \u' = u\ (3.24) 
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where now u 1 = du/dX. 

This scaling of x places us in a boundary layer next to the boundary at x = 1. 
Since the object of the rescaling is to bring back the highest derivative term at leading 
order, we see that we must choose 

5 = £, (3.25) 

and with this the equation can be written 

u" + v! = eu 2 . (3.26) 

The issue of appropriate boundary conditions for (3.26) is crucial to the procedure. 
One is fairly clear, that we specify 

u = 1 at X = 0; (3.27) 

this is just the right hand boundary condition. 

The other condition is less clear. We already have an outer expansion (3.17) which 
satisfies the left hand boundary condition, so the other condition can not be that. 
Nevertheless, we press on, aiming to find an inner expansion for (3.26) of the form 

u ~ U (X) + eU^X) + . . . , (3.28) 

and then, equating powers of e, we have 

uS + uL = o, 

U'{ + U[ = U 2 , (3.29) 

and so on. 

The leading order solution for Uq is 

U = e~ x + A (l - e~ x ) , (3.30) 

where A is a constant, and further constants arise in the determination of higher order 
terms. The issue is how we choose the constant A (and the other constants which 
arise at higher order). The idea is this: we have two expansions, an inner and outer 
expansion, each representing an approximation to the solution in different regions; one 
where 1 — x ~ 0(1), the other where 1 — x ~ 0(e). We suppose that each expansion 
should blend into the other in a fuzzy intermediate zone where e <C 1 — x <C 1. 
This region is called a matching region, and the method is thus called the method of 
matched asymptotic expansions. 

There is a methodical procedure for matching the expansions, which is beyond 
the scope of the present discussion, but is necessary if matching is to be carried out 
for the higher order terms. However, its implementation for the leading order terms 
is usually fairly straightforward. The simplest situation is when the outer solution 
tends to a constant at the boundary adjoining the boundary layer, while the inner 
solution tends to a constant as the inner variable becomes large. This is exactly the 
situation that applies here. In this case the matching condition is to ensure that 
these two constants are the same. Comparing (3.22) and (3.30), we see that the two 
solutions match if A = 0. Thus we have 

u~0, l-z~0(l), 
u~exp[-(l-x)/e], l-i~0(e). (3.31) 
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Uniform solutions 

Although the inner and outer expansions formally blend together, actual graphical 
illustration of them generally indicates discontinuity, and it is attractive to write a 
solution which is uniformly valid. It is possible to do this. The trick is to realise 
that, since both inner and outer expansions are supposed to be valid in the matching 
region, we can get a uniform solution by adding the inner and outer solutions and then 
subtracting that part of the solution which each represents in the matching region. 
In the present example this is trivial, since the common part of the outer and inner 
solutions is zero, and in fact the inner solution is uniformly valid. 1 

U 

1 

0.5 





0.5 1 1.5 

t 

Figure 3.1: Inner, outer and uniform solutions to (3.33) when e = 0.1. 

An easy illustration of this idea is the function 

u (t) = e -* - e ~ t/£ , (3.32) 
which is a uniform first order expansion of the solution to the initial value problem 

eii + u + u = 0, m(0) = 0, m(0) = -, e«l. (3.33) 

£ 

The first order outer solution of this problem is u = Ae~ l , while the first order inner 
solution is it = 1 — e~ T , where t = sT. Matching requires A = 1, and the common 
part of both approximations in the intermediate matching region is then just their 
common limit, which is one. Thus the uniform first order solution adds the inner 
and outer solutions and subtracts one. Figure 3.1 shows that both inner and outer 
approximations are well approximated in their regions of validity by the uniform 
expansion. 

J This lends substance to our earlier comment concerning exponential asymptotics. 
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3.2.1 Internal layers 



The question of which boundary or boundaries one should attach boundary layers to 
is a matter of possibilities. Because the rapidly varying solution in a boundary layer 
must blend into a slowly varying outer solution, it needs to decay in the matching 
region. Very often, this decay is exponential, and the choice of exponential decay 
rather than exponential growth determines the boundary layer location. For example, 
we found that (3.16) allowed a boundary layer at a right hand boundary, but (3.33) 
will only allow boundary layers on the left. 
More generally, the equation 

ey" + a(x)y' + b(x)y = (3.34) 

will allow left hand boundary layers where a > and right hand boundary layers 
where a < (see also question 3.6). This raises a potential problem, which we can 
illustrate by consideration of the nonlinear problem 

ey" + 2yy' + y(l-y 2 ) = 0, (3.35) 

with < £ « 1 and 

y(-l) = -l, y(l) = 1. (3.36) 

As we expect both the derivative terms to be important in any boundary layer 
(why?), we see that the sign of y at the boundaries precludes boundary layers there, 
and any region of rapid change must occur in the interior of the domain. Such 
regions are called interior layers, or in partial differential equation models they are 
called shocks (there will be more on these in chapter 7). 

(3.35) has outer solutions satisfying both boundary conditions given by 

y — ±1 for x ^ x , (3.37) 

where xq is indeterminate. There is also a possible outer solution y = 0, which we 
will ignore in this discussion. If the solution is unique, then it is odd (since — y(— x) is 
also a solution), and then x = 0. But because the problem is nonlinear, uniqueness 
(or even existence) is not guaranteed. 
In the transition region, we put 

x = Xq + eX- (3.38) 

then to leading order y satisfies 

y" + 2yy' = 0, y(±oo) = ±1, (3.39) 

and thus 

y = tanhX, (3.40) 

where a constant of integration which fixes where y = can be absorbed into the 
definition of xq. 



33 



This solution describes the interior layer, but there is nothing to determine the 
position xq of the transition. The reason for this is that the outer solutions y = ±1 are, 
like that of (3.16), accurate to all orders of e. Consequently, as far as the asymptotic 
solution is concerned, the location of the boundaries is immaterial, and it is because 
of this that the interior position of the transition layer is indeterminate. Its position 
(if a solution exists) presumably depends on the exponentially small correction terms 
which lie beyond all orders of the expansion in powers of e. 

A less degenerate case occurs if we modify the boundary conditions to be 

y(-l) = -B_, y(l) = B + , (3.41) 

where B± > 0. Suppose that B± < 1, and define 

5±=tanh,4±, (3.42) 

where A± > 0. The previous boundary condition (3.36) corresponds to the double 
limit A± — > oo. The outer solutions are y = or 

„ ± = ±tanh (I±£±^±). (3.43) 

Now at an internal transition layer, (3.39) allows solutions 

y = y*tanhy*X; (3.44) 

the jump is thus from — y* to y* , and this allows us to determine xq from (3.43): 

x = A + -A_. (3.45) 

Thus we obtain a self-consistent asymptotic solution, but only for a restricted range 
of boundary conditions. We can also see the degeneracy as A± — > oo. 

This example, apart from illustrating the occurrence of an interior layer, also 
shows how asymptotic solutions are not always straightforward, and often require 
guile to discover their intricacies. This is particularly true of nonlinear problems. 
To some extent, the problem arises here because the boundary conditions (3.36) are 
exact solutions of the equation. But it is also emblematic of the artificial nature of 
the example. It is more generally the case that good mathematical models of real 
physical systems do not have such pathological behaviour. 



3.3 Multiple scales and averaging 

A different kind of singular limit can occur if the domain of integration is infinite. This 
is typically the case in time dependent problems, and the usual sorts of illustrative 
example are models of oscillators. As an example, suppose that u(t) satisfies the 
initial value problem 

ii + u = £u 2 , u(0) = 1, «(0) = 0, (3.46) 
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where £ < 1, and u = du/dt. This is the equation of a weakly nonlinear oscillator, 
a spring with a weak amplitude-dependent spring constant oc 1 — eu. We seek a 
perturbation expansion of the form 

u ~ Mo + sui + . . . , (3.47) 

and substituting this into (3.46) we have the sequence of equations 

Mo + «o = 0, 
iii+ui = ul, 

U2 + U2 = 2u ui, (3.48) 

and so on. The initial conditions at t = for these are 

u = 1, u = 0, 
ui = 0, iii = 0, 

u 2 = 0, ii 2 = 0. (3.49) 

It is straightforward but tedious to solve these, and we find the resultant expansion 
for u to be 

u ~ cos t + £ \ — \ cos 2t — I cos t 

+£ 2 [-| + cos t + I cos 2t + i cos 3t + sint .... (3.50) 

This is a perfectly good asymptotic solution for bounded times, but the series becomes 
disordered at large times; specifically the third term becomes as large as the first 
when t = 0(l/e 2 ). This contradicts the assumption that the terms are decreasing 
asymptotically. 

The breakdown in the expansion is associated with the occurrence of the secular 
term tsint at 0(e 2 ), which arises through a resonance between the basic linear os- 
cillator and the forcing terms in the expansion. Now we know that the solutions are 
indeed oscillatory, and so we know that the amplitude cannot grow in this way. It is 
also clear why the method must break down. The period of the nonlinear oscillator is 
not constant, and depends weakly on the amplitude. But the unperturbed oscillator 
has fixed period, and our expansion method forces this period on the perturbation 
solution. It is this incorrect implicit assumption of an unperturbed period which 
generates the secular terms. 

This idea is simply illustrated. The function cos[(l — e)t] has period 2ir + 0(e). 
If we expand in powers of e, we find 

cos[(l — e)t] = cost + etsint + . . . ; (3.51) 

the secular terms arise through a misguided expansion of the oscillatory term. 
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3.3.1 Poincare-Lindstedt method 

There are various ways to fix this breakdown, of which the simplest is the Poincare- 
Lindstedt method. Essentially we need to generalise our naive perturbation method 
to allow for a drift in the period of the underlying oscillator. 

We do this by straining the time variable; we define a new time r via 

r = (1 + ewi + e 2 co 2 + . . .)*, (3-52) 

so that (3.46) becomes 

1 + 2eco 1 + e 2 (uf + 2co 2 ) + ...]«" + « = £« 2 , (3.53) 

where now u' = du/dr. The constants u>i are to be chosen precisely to eliminate the 
secular terms in the expansion. 

Expanding u as a perturbation series as before, we find 

u'q +u = 0, 

u'{ + Ui = Uq - W\Uq, 

u 2 + u 2 = 2uqU\ — {uj\ + 2oj 2 ^) u'q, (3.54) 

and we solve these successively as before. There are two little tricks which make it 
less algebraically nightmarish than our earlier calculation. One is to write the general 
solution of (3.54)i as 

u = Ae u + (cc), (3.55) 

where (cc) denotes the complex conjugate, and A is in general a complex amplitude. 
The second is not to specify A = |, but to leave it as an undetermined constant at this 
stage. The point is that we can in principle subsume all the fundamental solutions 
e lt which occur at each order into the definition of A, or equivalently we allow A to 
have its own expansion in powers of e. 
The 0(e) terms imply 

u'[ + Ul = A 2 e 2it + \A\ 2 + u x Ae {t + (cc); (3.56) 

terms oc e lt (or its complex conjugate) are secular producing terms, and so we must 
choose 

wi = 0; (3.57) 
there is no secular term at 0(e), and the (particular) solution is 

mi = \A\ 2 - lA 2 e 2it + (cc). (3.58) 

At the next order, we find (taking care to ensure we have all the nonlinear terms) 2 

u» + U2 = 2 [l\A\ 2 Ae« - lA 3 e 3it + co 2 Ae u + (cc)] ; (3.59) 



2 Note that (a + a)(b + b) — (a + a)b + (cc), where the overbar denotes the complex conjugate. 
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suppression of secular terms requires us to choose 

0J2 = -\\A\\ (3.60) 

and the (particular) solution is 

u 2 = ±A 3 e 3it + (cc). (3.61) 

As we surmised, the period P (in t) is weakly amplitude dependent, 

P = 2tt (l + \e 2 1 A\ 2 + ...), (3.62) 

and the coefficient ^ if \A\ = \ is consistent with the secular term in (3.50). This 
procedure can be carried on sequentially, but usually it is only the leading order 
frequency correction which is of interest, which is also why it is not really necessary 
to determine the amplitude of the fundamental beyond leading order. 



3.3.2 Van der Pol oscillator 

The Van der Pol oscillator is the equation for x(t), 

x + e (x 2 - l) x + x = 0, (3.63) 

and we suppose e is small and positive. Unlike the conservative nonlinear oscillator 
of the preceding section, there is a unique stable limit cycle of the Van der Pol 
equation. The Poincare-Lindstedt method will find this limit cycle, both its leading 
order amplitude and the frequency correction (see question 3.7), but to prove stability 
of the limit cycle, a more subtle approach is necessary. 

This retains the idea that the leading order solution will be of the form 

iwXe fl + (cc), (3.64) 

but rather than locking the time scale to a fixed modification of the original one, we 
instead allow the amplitude A to vary slowly in time. This is a generalisation of the 
Poincare-Lindstedt method, insofar as we can replace ^4exp[i(l + euj\ + . . .)t] (with 
constant A) by A(i)e lt , where the slow time t here could be taken to be et. Yet more 
generally, one can choose both a new strained fast time t* — t(l + . . .) as well as a 
slow time, and for strongly nonlinear equations (unlike the weakly nonlinear (3.63)), 
fast and slow times may have complicated functional dependence on t. At higher 
order, yet more slow time scales may be necessary. 

However, in its simplest form, such complications are not necessary. Consideration 
of (3.63) suggests that we define fast and slow times 

t* = cot, t = et, (3.65) 

where 

co = (l + e 2 co 2 + ■■■), (3.66) 
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and we then formally seek the solution as a function x(t*,t) of the two times t* and 
t. These give the method of multiple scales its name. 

The chain rule relates the t-derivatives to the partial t*- and t-derivatives, so that 
the equation (3.63) is embedded in the hyperbolic partial differential equation 



9 d 2 x d 2 x d 2 x / 9 \ / dx dx\ . „, 

w «a + 2e "W5i + £ V + e ( x ~ {"m- + s ai) + x = a (3 ' 67) 

We actually construct a perturbation solution of this equation in the two-dimensional 
(t*,t) space, but are only concerned with its value on the characteristic through the 
origin. There is some flexibility about how we impose initial conditions. If x(0) = xo 
and x(0) = vo, then the most natural initial data for (3.67) is to specify 

x = Xq on t = 0, 

Ox ~ 
oj-^— = v on t = 0. (3.68) 
at* 

As in our discussion of the Poincare-Lindstedt method, we shall not be much con- 
cerned with initial conditions. 

We now expand the solution in the usual way, 

x ~ x + ex\ + . . . , (3.69) 

where each Xi = Xi(t*,i). We then gain the following sequence of problems: 

d 2 x 

d 2 xi d 2 x / 2 x dx 

dt^ +Xl = ^Wdr^-^W (3J0) 

and luckily for us we need go no farther, since secular terms appear already at 0(e). 
(If they did not, then there would be no need of the 0(1/ e) slow time scale, and a 
slower time scale would be relevant.) The solution at leading order is 

x = A(t)e u ' + (cc), (3.71) 

and therefore the equation at the next order is 



d 2 Xl 
dt* 2 



+ x 1 = -2iA'e it * - iA 6 e 6ir + (2\A\ 2 - l)iAe lt ' - i\A{ 2 Ae lt ' + (cc). (3.72) 



3„3W i /r>l/l|2 iV- A At* ;\A\1aAV 



We do not even need to solve this, we simply need to choose A in order to remove 
the secular terms oc e lt * . We do this by requiring A to satisfy the Ginzburg-Landau 
equation 

2^ =A-\A\ 2 A. (3.73) 

(Jj v 

Multiplying this by A and its conjugate by A and adding yields the equation for \A\, 

= \A\ 2 - \A\\ (3.74) 
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which shows that \A\ — > 1 as t — > oo, and there is a stable limit cycle 

x~2cost. (3.75) 

The amplitude of the limit cycle is two, and the frequency is shifted by 0(e 2 ), but 
to find the magnitude of the frequency shift would require evaluation of the 0(e 2 ) 
terms. 



3.3.3 Averaging 

The method of averaging is another formal method which is particularly useful for 
weakly nonlinear oscillations, and also can be used rapidly to give a leading order 
evolution equation for the amplitude of motion. It relies on the fact that the unper- 
turbed system has a constant of integration, which will usually be the energy. We 
illustrate its use with the Van der Pol equation, 

x + e(x 2 - l)x + x = 0. (3.76) 

The unperturbed equation is that of simple harmonic motion, with constant amplitude 
a and phase velocity 9. This suggests the exact change of variables 

x + ix = ae id , (3.77) 

as a consequence of which 

The energy |a 2 is conserved when e = 0, and therefore we can expect that a will be 
slowly varying in the motion. Specifically, we have from (3.76) 

ad = -e{x 2 - l)x 2 ; (3.79) 

using (3.77), this implies 

a = ea(l - a 2 cos 2 0) sin 2 9, (3.80) 

and we can also deduce that 

9 = -1 + 0(e), (3.81) 

though this is of less concern. The central idea of the averaging method is that 
since a and 9 are approximately over the underlying period of the oscillation, we can 
effectively integrate the amplitude equation over this period keeping a constant. The 
result of doing this to (3.80) is to take the average (in 9) of the right hand side, and 
this yields the approximate amplitude equation 

a w ea (| - |a 2 ) . (3.82) 

Putting a = 2A, we have the same evolution equation as (3.73), and 

a 2 as t oo, (3.83) 

as we found before. 
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3.4 Notes and references 



There is a wealth of books on perturbation methods. A rapid assessment is that by 
Hinch (1991); longer classics are those by Kevorkian and Cole (1981) and Bender and 
Orszag (1978); the last of these annotates its text with many illustrative numerical 
examples. 

Exercises 

3.1 Find two terms of a regular perturbation expansion of the form 

oo 


to the boundary value problem 

u" + u = eu 2 , u(0) = 0, «(7r/2) = 1, 
when e is small and positive. 

Show by induction that u n can be written in the form 

n+l 

U n = J2 a k eikX - 
-(n+l) 

Repeat the exercise for the equation 

u" + u = eu 3 , 

with the same boundary conditions. Can you modify the inductive step to find 
a general form of expression for «„? 

3.2 Suppose that x(t) and y(t) are real and satisfy the pair of first order equations 

x = —y — 2eaxy, 
y = x + e{x 2 -y 2 ), 

in which a is constant, and x and y satisfy initial conditions 

x(0) = l, 2/(0) = 0. 

By eliminating y, find a second order equation for x of the form 

x + x = ef(x, x; a), 

and give an expression for the function /. Write down the associated initial 
conditions for x. 

Use perturbation theory to find two terms of an expansion for x in powers of e. 

Find an exact solution for the case a = 1 by writing z = x + iy, and show 
that the expansion of this solution in powers of e agrees with the perturbation 
expansion. Deduce that in this case the perturbation expansion converges for 
sufficiently small e. 
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3.3 The functions u(t) and v(t) satisfy the initial value problem 



eu = v — u, 

v = —u, 

with the initial conditions 

u(0) = 0, v(0) = 1. 

Show that if the outer solution is chosen to satisfy the initial condition for v, 
then this outer solution is v = e~ l . Show that there must be a boundary layer 
near t — 0, and by rescaling t, write down a suitable version of the equations in 
this boundary layer, and find a leading order approximation for u and v in the 
boundary layer. 

Hence write down a uniform approximation for u. 

3.4 The Michaelis-Menten model of enzyme kinetics is the pair of equations 

8 = -s(l - c) + (K - L)c, 
EC = s — (s + K)c, 

where e <C 1 and K,L = 0(1). The initial conditions are that 

s(0) = 1, c(0) = 0, 

and the constants e, K and L are all positive. Find leading order approximations 
to the solution. 

3.5 Find leading order approximations (including boundary layers, where necessary) 
to the solutions of the boundary value problems 

ey" + 2y' + e y = 0, y(0) = y(l) = 0, 
ey" - x 2 y' - y = 0, y(0) = y(l) = 1. 

3.6 Suppose that y(x) satisfies the boundary value problem 

ey" + a(x)y' + b(x)y = 0, 

with < e < 1 and 

2/(0) = A, y(l) = B. 

Suppose also that a(x) and b(x) are analytic in [0,1] (i.e., have convergent 
Taylor series about any point in this interval). 

(i) If a > 0, find an approximate solution and show that it has a boundary 
layer at x = 0. 
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(ii) If a < 0, find an approximate solution and show that it has a boundary 
layer at x = 1. 

(iii) Finally, if a < for x < x and a > for x > x , where < x < 1, 
show that no boundary layers at the end points can exist, and therefore that 
an interior layer must exist at xq. 

Suppose that 

b(x ) 
P a'(x o y 

Show that as x — > xo±, the outer solutions in x ^ xo satisfy 

y ~ c ± |x - xol^, 

where the constants c± should be specified. 
Hence show that by rescaling x and y as 

y = { £ a'(x )f /2 y, x = x + {ea\x )} l l 2 'X, 

the equation can be approximately written in the transition region as 

Y" + XY' + @Y = 0, 

and explain why suitable matching conditions for the solution of this equation 
are 

Y ~ c±\X\ p as X ±oo. 

3.7 The Van der Pol oscillator is represented by the equation 

x + £ (x 2 - X + X = 0, 
and we assume that < e -C 1. 

Use the Poincare-Lindstedt method to find the amplitude and frequency cor- 
rection of the periodic solution of the equation. 

3.8 Duffing's equation is 

x + x + ex 3 = 0. 

By first finding a first integral, show that all solutions are periodic if e > 0. If 
£ <C 1, use the Poincare-Lindstedt method to construct approximate solutions, 
and find the frequency correction of a periodic solution in terms of e and its 
amplitude A. 
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Chapter 4 

Stability and oscillations 



If we move from first order systems to second order systems of the form 



x 



y 



f(x,y), 
g(x,y), 



(4.1) 



more interesting phenomena can occur. There is indeed a hierarchy of complexity 
which one ascends as the order of the equation increases. As we saw in the preceding 
chapter, first order equations have steady state solutions which are alternately stable 
and unstable, and the instability is direct, in the sense that loss of stability as a 
parameter changes leads to a transient migration towards another fixed point. 

In second order systems, a different kind of instability can occur. As well as the 
direct instability, or exchange of stability, between different fixed points, oscillatory 
instability can occur, and the consequence of such instabilities is that permanent oscil- 
latory (periodic) solutions may exist: in two dimensions there is dynamic behaviour. 1 

4.1 Linear stability 

We illustrate the technique of linear stability analysis for (4.1). The analysis applies 
(and is easy) in two dimensions, but evidently the method applies in n dimensions, 
although the classification of behaviour takes its essence from the two-dimensional 
example. 

A steady state of (4.1) is a constant pair (x ,J/o) which satisfies the equations, 
i.e., f(xo,yo) = g(xo,yo) = 0. For small perturbations about this state, we write 



x = x + X, y = y + Y, 



(4.2) 



and linearise the system (4.1) to obtain the approximate equations 




(4.3) 



1 When one proceeds to third or higher order systems, more exotic behaviour can occur: period- 
doubling, quasi-periodic oscillations, and chaos, for example. 
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where the partial derivatives are evaluated at (x ,y ). The matrix 

M = ( fx f y \ / 44 n 

\9x 9y J 

is sometimes called the community matrix (particularly in applications in population 
biology), and its trace and determinant will be denoted 

T = tr M, D = det M. (4.5) 

Because M is a constant matrix, solutions of (4.3) are of the form 

; A * =ue xt , (4.6) 



Y 



where u is an eigenvector of M and A is the corresponding eigenvalue; thus A is a 
root of the quadratic equation 

A 2 - TA + D = 0. (4.7) 

Classification of the different kinds of behaviour follows from the different combina- 
tions of pairs of values of A. For D > T 2 /4, the roots are complex, and the consequent 
phase portrait is a spiral; unstable if T > 0, and stable if T < 0. Examples are shown 
in figure 4.1. 



Figure 4.1: Unstable (left) and stable (right) spirals solving (4.3) with 
M = ( _\ _q 4 ) (left) and M = ^ ^ ) (right). 

If < D < T 2 /4, the fixed point is a node, with two real eigenvalues of M having 
the same sign. The node is unstable if T > and stable if T < 0. An example is 
shown in figure 4.2. 

Finally, a saddle point is shown in figure 4.3. This occurs if D < 0. The eigenvalues 
are real and opposite in sign. There is one stable direction and one unstable, and 
the fixed point is thus unstable. Putting these results together, we see that stability 
occurs if and only if D > and T < 0, and the classification of fixed points as spiral, 
node or saddle can be represented in (D, T) space as shown in figure 4.4. 
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Figure 4.2: Unstable node with M 
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Figure 4.4: T-D parameter space indicating location of stable and unstable spirals 
and nodes, and saddles. 

4.2 Linear systems 

It is easy to extend the discussion of two-dimensional systems to systems of n dimen- 
sions. We consider first autonomous systems, i. e., of the form 

x = Ax, (4.8) 

where x G R" and A is an n x n constant matrix. 2 Problems of this type arise when 
studying the stability of steady states of n-th order ordinary differential equations of 
the form x = fix), when A = Df(xo) is the Jacobian of / at the fixed poinbt xq. As 
before, solutions exist of the form 

x = e xt u, Au = Am, (4.9) 

and solutions grow or decay according as Re A ^ 0. x = is stable if Re A < for all 
eigenvalues A, but is unstable if at least one A has positive real part. 

A succinct representation of the general solution follows from defining the expo- 
nential matrix e tA . This slightly strange looking beast is defined in the obvious way, 
as 

t 2 A 2 

e tA = I + tA+ — + .... (4.10) 

For this definition to make sense, we need to borrow some ideas from linear algebra 
and functional analysis. We need to define convergence in the space ofnxn matrices, 
which requires the definition of an appropriate norm. Matrices being overgrown 

2 We are straying into the territory of dynamical systems, in which it is conventional not to use 
the bold x for vectors, and we will follow this convention in this section. 
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vectors, these are supplied in the usual way. Convergence of the exponential series is 
then automatic, and indeed the usual exponential laws apply, in particular e M e N = 

e M+N for matrices M and N, and thus also (e M ) = e~ tA . 
The general solution of (4.8) can be written 

x = e tA x , (4.11) 

where x = xq at t = 0. That this is a solution follows from first principles, since 

the definition of e tA implies that — e tA = Ae tA , and Picard's theorem implies it is the 

at 

only solution. 

With the exponential matrix in hand, it is easy to find the solution of the inho- 
mogeneous problem 

x = Ax + g(t). (4.12) 



Noting that 

it quickly follows that 



d [e~ tA x] = -e~ tA Ax + e~ tA x, (4.13) 



dt 



x = e tA x + 



f t e^ A g(s)ds. (4.14) 
Jo 



4.2.1 The fundamental solution 

We now generalise the preceding discussion to non-autonomous systems of the form 
(4.8), in which A = A(t). Let x^ denote the solution of this with 

*?'«»=*> = { J; i;t, (^5) 

5ik is the Kronecker delta. The vectors x^ are initially linearly independent, and 
the matrix whose columns are these vectors is called the fundamental matrix, or 
fundamental solution, $(£), and $^ = x\ k \ In particular, 

$ = A$, $(0) = /. (4.16) 

The fundamental matrix plays the role of (and, indeed, is) the exponential matrix 
of an autonomous system. It is easy to show that the solution of the inhomogeneous 
equation (4.12) is 

$(t)x + [ ^(m- 1 (s)g(s)ds. (4.17) 



X 



4.2.2 Floquet theory 

Floquet theory concerns the case where the matrix A(t) is periodic. Without loss of 
generality we will suppose that its period is 2ir. Such systems arise when considering 
the stability of periodic orbits. Suppose, for example, that the system of ordinary 



47 



differential equations x = f(x) has a 2-7r-periodic solution x p (t). To study its stability, 
we put x = x p + X , and linearise the equation for small X. We obtain 

X = A(t)X, (4.18) 

where A = Df{x p (t)} is the Jacobian of / on the periodic orbit, and is evidently 
periodic. The question we ask is whether the solutions X grow or decay in time. 

It is not possible to provide such a good recipe as when A is constant (calculate the 
eigenvalues of A), but thanks to Floquet's theorem we can come close. Again, we let 
$ be the fundamental solution for A, thus (4.16) applies. We define the monodromy 
operator to be the matrix 

M = $(2tt); (4.19) 

in the language of dynamical systems it is the Poincare map on the surface of section 
t = mod27r. The difficulty with Floquet theory lies precisely in determining this 
map. As we shall see, stability is directly determined by the eigenvalues of M. 

We need to show that zero is not an eigenvalue of M. To do this, we define the 
Wronskian W of the system to be the determinant of From first principles, we can 
deduce Liouville's theorem, that 

W = (tvA)W, (4.20) 

and it follows from this that W is never zero (since W = 1 at t = 0). This implies, 
in particular, that det M ^ 0, and since the determinant of a matrix is the product 
of its eigenvalues, we see that the eigenvalues of M are non-zero. 

If the eigenvalues of M are /j, s , we can therefore define quantities A s by 

Hs = e 2 * x % (4.21) 

and we note that Re A > if and only if \p\ > 1. Next we define the diagonal matrix 

A = diag(A s ), (4.22) 

and we define a matrix \J>(t) of fundamental type by 

ijr = A*, *(0) = /, (4.23) 

Now suppose that M is diagonalisable, as is normally the case, thus there is a 
matrix P such that 

P~ X MP = diag (e 2nXs ) . (4.24) 

We then have that 

Ptf (0) = P = $(0)P, (4.25) 

and 

P$(2tt) = $(2tt)P, (4.26) 
and it follows from this that the matrix 

B{t) = §P^- 1 P- 1 (4.27) 
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is 27r-periodic. Further, if we define X = B(t)Y, then 

Y = (PAP" 1 ) Y. (4.28) 

Therefore the stability of the periodic orbit x = x p (t) is determined by the sign of 
the real part of the eigenvalues X s of A, or equivalently by whether the characteristic 
multipliers /j, s of M lie inside or outside the unit circle in the complex plane. Unfor- 
tunately, these must in general be computed numerically, since M must be computed 
numerically. On the other hand, this is a relatively simple computational task, if 
the periodic orbit can be computed. Of course, the computation itself then tells us 
directly whether the periodic orbit is stable or not. 



4.2.3 Mathieu equation 

4.3 Nonlinear stability 

4.4 Phase plane analysis 

To go beyond linear stability analysis and complete the phase portrait of solution 
trajectories is the subject of phase plane analysis, and the most interesting feature 
of the phase plane is that periodic oscillations can occur. An illuminating example is 
illustrated in figure 4.5, and is typified by (but is not restricted to) the equations 

x = g{x) - y, 

y = y-h(x), (4.29) 

where the functions g and h are as shown in the figure: g is unimodal (e.g., like 
g = xe~ x ) and h is monotonic decreasing (e.g., like h = l/(x — c)). The graphs of 
g(x) and h(x) (and more generally, the curves where x = and y = 0) are called 
the nullclines of x and y, and it is simple to see that where they intersect, there is 
a steady state solution, and also that in the four regions separated by the nullclines, 
the trajectories wind round the fixed point in an anti-clockwise manner. 

The next issue is whether the fixed point is unstable. If we denote it as (x*, y*), 
write x = x* + X, y = y* + Y, and linearise for small X and Y, then 

(4.30) 

where the derivatives are evaluated at the fixed point. The stability of this two- 
by-two system with community matrix A = ^ ^ ^ is governed by the trace 

and determinant of A, as indicated in figure 4.4. In the present case, tr A = g' + 1, 
det^4 = g' — h 1 , so that for the situation shown in figure 4.5, where h! < g' < 0, 
del A > 0, and the fixed point is an unstable spiral (or node) if g' > —1. When 
g' = —1, there is a Hopf bifurcation, and if the system has bounded trajectories (as is 
normal for a model of a physical process) then one expects a stable periodic solution 
to exist. Figure 4.6 illustrates an example. 




49 




Figure 4.5: Nullclines for (4.29). 




Figure 4.6: Typical form of limit cycle for a system with nullclines as in figure 4.5. 
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4.4.1 An example of Murray 



Murray (2002, p. 116) gives an example of a population interaction model repre- 
senting aquatic populations of herbivores (e.g., zooplankton) and (phyto-)plankton. 
Denoting the populations as x (phytoplankton) and y (herbivore), the model is 



x 



y 



x 



k — x — 



y 



ay 



x 



1 + x 



1 + x 
ay 



(4.31) 



and as normal in population models we are concerned only with the first quadrant 
x > 0, y > 0. We suppose the constants k, a and a are positive. 




Figure 4.7: x and y nullclines for (4.31), with a = 1 and k = |. Arrows indicate the 
direction of the trajectories in the four sectors bounded by the nullclines. 

Define the functions 

F(x) = (k-x)(l + x), G(x) = ^-^. (4.32) 

The x and y nullclines are y = F(x) and y = G(x)/a, respectively, and a typical 
situation is shown in figure 4.7. There are two fixed points at (0,0) and (k,0), and 
there is always a fixed point (xo, yo) in the positive quadrant. If k < 1 as in the figure, 
then this fixed point is uniquely defined as there is only one intersection point of the 
nullclines (in which both x and y are positive). 

The stability at the origin can be found by expanding for small x and y, whence 
we find 

x & kx, y ~ ay(x — ay). (4.33) 

Note that the origin is degenerate in the sense that the linear approximation for y is 
just y ~ 0. It is clear from the behaviour of the solutions of (4.33) that the origin is 
a kind of saddle point. It is similarly easy to show that (k, 0) is a saddle. 
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To find the stability of the other fixed point, we write 

X = X + X, y = y +Y, (4.34) 

and linearise, so that 

(?H(?)< (436) 

where we use the fact that at (xo,yo), ayo = xq/(1 + xo); F' and G' are evaluated at 

Xq. 





Figure 4.8: Phase plane trajectories for (4.31), with a = 1, a = 0.4 and k = 0.5. 
From this we find that 

trM = ayo(F'-a), 
detM = aayl[G' - aF'}. (4.36) 

It is evident from the way the nullclines intersect (since their slopes are just G' /a 
and F') that for the situation shown in figure 4.7, det M > 0. Therefore stability is 
associated with the sign of trM, and in particular there is a Hopf bifurcation to a 
periodic orbit if F' > a. For the situation of figure 4.7, where k < 1 and thus F' < 0, 
we see that this fixed point is a stable node or spiral. Figure 4.8 shows this behaviour. 

The possibility of instability and the consequent occurrence of a periodic oscilla- 
tion occurs if F'{xq) > 0. Since xq depends only on a and k, we see that in this case, 

there will be instability if a is small enough. The question then arises, for what a 

k — 1 

and k is F'(xo) > 0? Since F'(xq) = k — 1 — 2xo, we need xq < — - — (remember we 

also must have k > 1). Since F(x ) = G(x )/a and G is monotonically increasing, 
we thus also have 

F{xo) = «W K I , (4 . 37) 

a a V k + 1 1 
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and since F(x ) < F(k — 1) = \(k + l) 2 , we see that a necessary condition for 
instablity to occur is that 

4(Jb - 1) 

0<a< (¥TTF < 438 > 

When k is close to 1, the upper boundary of this region becomes a good estimate 
of the true upper limit. Figure 4.9 shows an example of periodic behaviour when 
a = 0.05, a = 0.9, and k = 3. 




Figure 4.9: A stable periodic solution when a = 0.05, a = 0.9, and k = 3. The 
nullclines are shown as the dashed lines. 
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Figure 4.10: Relaxational solution of (4.31) when a = 0.02, a = 0.9, and k = 3. The 
graph is that of x(t). 



It is evident from this phase portrait that the limit cycle has a relaxational char- 
acter, in the sense that there is a (longish) period when x is very small. Of the 
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parameters used in figure 4.9, we see that a is small, suggesting that the limit a — > 
is responsible for this behaviour. This is confirmed in figure 4.10, which shows the 
time sequence of x(t) when a = 0.02, a = 0.9, and k = 3. It is not uncommon to see 
near extinction oscillations in population models. It is possible, though not entirely 
straightforward, to analyse such oscillations using singular perturbation methods. 

Multiple steady states 

There is always at least one steady state with y > 0, but inspection of F(x) and G(x) 
shows that there can be more than one. Figure 4.11 shows an example. Preceding 
discussion indicates that if we label the three equilibria in this case as (xj,yj), with 
i = 0,1,2 and x < %\ < %2, then the middle steady state (xi,yi) is always an unstable 
saddle (detM < 0), the upper steady state (2:2,2/2) * s always stable (detM > and 
trM < 0), and the stability of (x , J/o) is determined as before. 




Figure 4.11: Nullclines of x and y when k = 10 and a = 0.027. 

It is clear from figure 4.11 that, since F increases with k, there will be an interval 
of fc-values, k_ < k < k + , for which multiple steady states occur. The values k± can 
be evaluated by the condition that there be a double root of F(x) = G(x)/a, and this 
leads to the parametric prescription for k± as 

2x 2 x-1 

k = —r a= (TTW (4 ' 39> 

which describes a narrow sector in the (k,a) plane, as shown in figure 4.12. The 
upper and lower boundaries of this sector are given, from (4.39), by 

1 4 

a ~ — , a ~ — r as k — > 00. (4.40) 
4k k 2 

Figure 4.12 suggests that oscillations do not occur when there are multiple steady 
states. 
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Figure 4.12: The regions in which instability of (x ,Z/o) ma y occur (U) and in which 
there are multiple steady states (M). The figure suggests (but does not prove) that 
when multiple steady states occur, both the upper and lower branches are stable. 

4.4.2 Travelling waves on the Martian north polar ice cap 

Phase plane analyses often arise when one seeks travelling wave solutions of partial 
differential equations, as will be discussed further in chapter 7. An example of this 
approach occurs in the study of a model describing accretion and sublimation of ice 
on the north polar ice cap of Mars. As can be seen in figure 4.13, a series of arcuate 
canyons occur on this ice cap, marked by dust deposits on the steep scarps. The ice 
cap is some three kilometres deep and 1000 kilometres wide, while the canyons are of 
typical depths 500 m, and spaced some 50 km apart. Their occurrence is mysterious, 
and unlike anything seen on the ice sheets of Earth. 

One model which has been studied with a view to explaining them monitors the 
effect of atmospheric dust concentration on sublimation rate. To explain the spiral 
wave like form of the troughs, travelling wave solutions of the model are sought, and 
these are described by 



in which q = q(c, s); s represents surface slope, q is the net rate of sublimation, and c 
is the atmospheric dust concentration. The prime denotes a spatial derivative in the 
moving frame of the travelling wave, v denotes the travelling wave speed, while u is 
the atmospheris wind speed, and <j> is the dust fraction of the ice, and these are all 
taken as constants; u and </> are positive, but v is unknown. 
The sublimation function q is taken to be 




= vs-q, 
= H, 



(4.41) 



Q = /(c) + ma, 



(4.42) 



55 




Figure 4.13: The spiral canyons of the Martian north polar ice cap. 

with m being positive and /(c) having a typical cubic shape exemplified by 

/( c ) = ( c -l)( c -2)(l-^). (4.43) 

The model has three steady states where s = 0, and c = 1, c = 2 or c = M; 
we denote these steady states as P, Q, R, respectively, as shown in figure 4.14. The 
results of a phase plane analysis are then the following. When < v < m, P and R 
are saddles, while Q is a node or spiral, which undergoes a Hopf bifurcation at 

1 £1 

v = Vc = m -ZL j (4.44) 

assuming v c is positive. Q is unstable for v > v c , and numerical computations suggest 
that the Hopf bifurcation is supercritical, so that there is a stable periodic solution 
for v > v c . This orbit exists up to a value v = Vh, at which it becomes heteroclinic to 
P, and disappears in a 'blue skies' bifurcation. The stable orbit and its disappearance 
are shown in figures 4.14 and 4.15. 

If v > m, P and R are still saddles, and Q is an unstable node or spiral. It is 
possible that the limit cycle may still exist, and this probably depends on the size of 
/; it is also possible for the limit cycle to disappear by coalescence with R, but this 
seems unlikely for large M. In fact, it is found that as M —> oo in (4.43), the limit 
cycle remains finite and independent of M. 

4.5 Relaxation oscillations 

It is a general precept of the applied mathematician that there are three kinds of 
numbers: small, large, and of order one. And the chances of a number being 0(1) are 
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Figure 4.14: Approach to stable limit cycle solution of (4.41) using (4.43) with M = 
1000, and v = 0.05, </> = 0.2, m = 0.14, and u = 2. The parabolic curve is the 
nullcline q = 0. 
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Figure 4.15: Disappearance of limit cycle in heteroclinic connection with P. Param- 
eters as for figure 4.15, except that v = 0.06. 
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Figure 4.16: Typical form of relaxation oscillation in phase plane for (4.45). 

not great. Thus for systems of the form (4.1), it is often the case in practice that the 
time scales for each equation are different, so that in suitable dimensionless units, a 
second order system might take the form 



where the parameter e is small. The system (4.45) is essentially the same as (4.29) 
with time reversed, but now suppose that the nullclines are as shown in figure 4.16, i.e. 
g has a cubic shape. Trajectories now rotate clockwise, and linearisation about the 
fixed point yields a community matrix A with tr A = — (g'/e) — 1, det A = (g 1 — h') /e, 
thus with g' > h', the fixed point is a spiral or node, and with e <C 1, tr A px —g'/e > 0, 
so it is unstable. Thus we expect a limit cycle, and because £ C 1, this takes the 
form of a relaxation oscillation in which the trajectory jumps rapidly backwards and 
forwards between branches of the x nullcline. For e <C 1, x rapidly jumps to its 
quasi-equilibrium y ph g(x), and then y migrates slowly (x [h(x) — g(x)]/g'(x)) 
until g' = and x jumps rapidly to the other branch of g. Figure 4.17 shows the time 
series of the resulting oscillation. The motion is called 'relaxational' because the fast 
variable x 'relaxes' rapidly to a quasi-stationary state after each transient excursion. 

4.5.1 The Van der Pol oscillator 

The relaxational form of the van der Pol oscillator is 



ex 



y-g(x), 

h(x) - y, 



(4.45) 



ex + (x 2 - l)x + x = 0, £«1. 



(4.46) 
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Figure 4.17: Time series for x corresponding to figure 4.16. 

This form of the equation is appropriate when the 'damping' (first derivative) term 
is large. If we define 

y = ex + |x 3 - x, (4.47) 
then the appropriate phase plane equations are 

ex = y — S(x), 

V = -x, (4.48) 

where 

S(x) = |x 3 - x. (4.49) 

The x-nullcline is the cubic y = S(x); above this, trajectories move rapidly to the 
right, and below it they move rapidly to the left. Thus, from an arbitrary initial 
location, a trajectory will rapidly approach either upward sloping part of the slow 
manifold y = S(x), on which the motion is slow. On the right side, in x > 1, y 
decreases until it reaches the base of the curve at x = 1, where it undergoes a rapid 
transition (an interior layer, in fact) to the left side of the curve in x < —1. There, 
it increases until a symmetric transition occurs at x = — 1. In this way, a periodic 
oscillation occurs. The typical form of the trajectories is shown in figure 4.18, which 
shows that the relaxational form is already appropriate when e = 0.1. The time series 
of x versus t is very similar to figure 4.17. 

It is possible to derive characteristics of the oscillation very simply from this 
description. The maximum amplitude is x = 2, following which there is a slow 
decline to x = 1. Because of the symmetry of the oscillation, the time to get from 
x = 2 to x = 1 is approximately P/2, where P is the period (since the fast transitions 
between are very short). Thus the period is 

P re 2 f =1 dt « f 2 dX « 3 - 2 In 2. (4.50) 

Jx=2 Jl X 
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Figure 4.18: Trajectories of (4.46) approaching the limit cycle (thick line) when e = 
0.1. Also indicated is the nullcline for x given by y = S(x). 

This follows because on the slow manifold, dt ?s —dy/x and y « S(x). 

4.6 Belousov-Zhabotinskii reaction 

4.7 Notes and references 

Phase plane methods are described in many books on ordinary differential equations, 
for example that of Jordan and Smith (1999). 

The blue skies catastrophe and various other bifurcations of exotic type are de- 
scribed in a number of books on dynamical systems, for example that of Thompson 
and Stewart (1986). 

The model of formation of the the spiral troughs on the Mars north polar ice cap 
is the subject of a dissertation by Zammett (2004). 

Linear systems of ordinary differential equations are the stuff of many basic text- 
books on the subject. A particularly clear and direct discussion is that by Arn'old 
(1973). 

Exercises 

4.1 « and v satisfy the ordinary differential equations 

u = k\ — + ksu 2 v, 
v = &4 — k 3 u 2 v, 
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where ki > 0. By suitably scaling the equations, show that these can be written 
in the dimensionless form 

u = a — u + u 2 v, 
v = b — u 2 v, 

where a and b should be defined. Show that if u, v are initially positive, they 
remain so. Draw the nullclines in the positive quadrant, show that there is a 
unique steady state and examine its stability. Are periodic solutions likely to 
exist? 

4.2 Suppose that 

X = A(t)X, (*) 

where X e R™ and A is 27r-perodic. Define what is meant by the fundamental 
matrix and write down the solution of (*) satisfying X = X at t = 0. 

Define what is meant by the monodromy matrix M, and show that it has no 
zero eigenvalues. Describe how it can be used to prove Floquet's theorem, that 
there is a 27r-periodic matrix B(t) such that X = BY, and 

Y = CY, 

where C is a constant matrix. Explain how the eigenvalues of C are related to 
those of M. 

4.3 Use a phase plane analysis (with y = x) to analyse the solutions of the nonlinear 
oscillator 

x + V'(x) = 0, 

where 

V(x) = \x 2 - \x A . 
In particular, show that the origin is genuinely a centre. 
Construct the phase planes for the systems 

(i) x = y, y = -V'(x) - ey, 

(ii) x = y, y = -V'(x) - ey 3 , 

where e is small. 
Show that if 

x = y, y = -V'(x)-sD{E)y, 

where 

E=\y 2 + V{x), D(E) = E-l 

then there is a stable limit cycle in the region < E < |, and draw the 
corresponding phase plane. 
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4.4 An energy balance model for the average atmospheric (absolute) temperature 
T of the Earth is given by 



pc p hf = \Q(1 - a) - a-fT 4 , 



where p is air density, c p is the specific heat, h is the effective height of the 
atmosphere, Q is the received short wave solar radiation, a is the albedo coef- 
ficient (reflectivity) of the Earth's surface, a is the Stefan-Boltzmann constant, 
and 7 < 1 is a numerical factor which represents the blanketing (greenhouse) 
effect of the atmosphere. 

Non-dimensionalise the equation using a temperature scale To and a time scale 
to, so that it takes the dimensionless form 



Show that if the observed albedo is 0.3 and the observed (steady state) temper- 
ature is 288 K, then T = 315 K. 

Use the values Q = 1370 W m" 2 and a = 5.67 x 10" 8 W m" 2 K" 4 to find what 
the greenhouse factor 7 must be to obtain this value of T . 

Use this value of 7, and the values p = 1 kg m -3 , c p = 10 3 J kg -1 K _1 , h = 10 4 
m, to compute the adjustment time scale to- 

The albedo of the Earth varies principally due to the extent of sea and land 
ice, which in itself depends on temperature. Sea ice varies seasonally, but land- 
based ice sheets grow over time scales of thousands of years. Suppose that as 
a consequence the albedo adjusts towards an equilibrium value A(T) on a time 
scale ti ~ 1000 y (years), and is modelled by the equation 



where v = to /it- 

We can expect A(T) to be a monotonically decreasing function, with, say, A w 
0.8 for very low T, and A w 0.2 for very high T. If A changes rapidly enough 
with T in an appropriate range, show that it is possible to have three steady 
states, and construct the phase plane for (T, a) in this case. Hence show that 
the smallest and largest steady states of T are stable, and the intermediate state 
is unstable. 

4.5 The sizes of two populations, one being a predator and the other its prey, are 
described by the equations 



f= 1-a-T 4 . 



Ua = A(T) - a. 



Show that the non-dimensional form of this equation is 



a = v[A(T) - a], 



u(l — u) 



auv 



u 



u + d' 
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where a, b and d are positive constants. 

Show that if u and v are non-negative, then they remain so. 

Show that there is a unique steady state in the first quadrant, and show that 
there is a region in (d, a, b) space where this steady state is unstable. 

Draw the possible form(s) of the first quadrant phase plane in this case. 

4.6 A nonlinear oscillator is given by the equation 

ex + (x 4 - l)x + x = 0, £«1. 

A suitable phase plane is spanned by (x, y), where y = ex + ^x 5 —x. Describe the 
motion in this phase plane, and find, approximately, the period of the relaxation 
oscillation. What happens if e < 0? 

4.7 The Belousov-Zhabotinskii chemical reaction can be approximately described 
by the two component pair of ordinary differential equations 

Z = -yX - Z, 

where e and 5 are very small, and 7 is 0(1). Show that relaxation oscillations 
will occur for 7 within a certain range (7., 7+), and give approximations for the 
values of 7±. 
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Chapter 5 



Oscillation and comparison 
methods 

The study of the stability of fixed points of systems of ordinary differential equations 
leads to the computation of the growth rate exponents a as eigenvalues of matrices. In 
a similar way, the computation of the frequencies of small oscillations of Hamiltonian 
systems leads to algebraic eigenvalue problems. 

Both of these situations arise because of the underlying linearity of the system, 
which allows the possibility of finding separable solutions. In partial differential equa- 
tions, the same situation arises, but now the computation of eigenvalues requires the 
solution of differential, rather than algebraic, equations, and the corresponding solu- 
tions of these are called eigenfunctions (as opposed to eigenvectors). In this chapter 
we discuss properties of these eigenvalue problems. 

5.1 Derivation of Sturm-Liouville systems 

We begin by giving some examples of eigenvalue problems which arise through the 
technique of separation of variables in partial differential equations. All of these 
examples take the form of (second order) Sturm-Liouville systems, that is, ordinary 
differential equations for y(x) of the form 

(py')' + (Xr - q)y = 0, (5.1) 

together with suitable homogeneous boundary conditions of the general form 

y' + ay = at x = a,b. (5.2) 

In this chapter we will mostly take a = oo, thus 

y = at x = a,b. (5.3) 

Sturm-Liouville systems possess the functional analytic property of being self- 
adjoint (the definition of which will appear in due course). Indeed, our discussion 
of Sturm-Liouville systems will frequently skirmish with concepts from functional 
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analysis, but the material will be presented more in the manner of a tour, and will 
not dwell overmuch on the niceties of detailed proofs. For these, reference must be 
made to any number of classical textbooks (see section 5.4). 

5.1.1 Waves in a bucket 

The first example we consider is that of free surface fluid flow in a cylindrical container: 
waves in a bucket, or a petri dish, or a pond. We suppose that the fluid, when 
undisturbed, occupies the region —h < z < 0, < r < a, where r, z (and 6) 
are cylindrical coordinates. When the fluid is disturbed, the upper (free) surface is 
perturbed to a surface z = rj(r,9,t), whose location forms part of the problem to be 
solved: it is therefore a free boundary problem. 

For irrotational flow, there is a velocity potential </>, which satisfies Laplace's 
equation 

V 2 = <j) rr + V + ^<j)ee + <t>zz = (5.4) 
together with the boundary conditions of no flow through the base and the sides: 

4> r = on r = a, 

(p z = on z = —h. (5.5) 

On the top surface, the kinematic and Bernoulli conditions imply 

<Pt + gv = 0, (p z = Vt on z = 0, (5.6) 

where it is assumed that the free surface elevation ij is small, so that the boundary 
conditions can be linearised about z = 0. In all these equations, subscripts indicate 
partial differentiation. 

The equation and boundary conditions are linear, and are autonomous in t, 9 and 
z (though not r), so that separable solutions can be found. These take the form 

= R( r )e iuit+ime cosh{k(z + h)}, rj = A(r)e iujt+ime , (5.7) 

where, because of the linearity, we allow R and A to be complex, and it is understood 
that we may take the real part of the solutions. Note that we require m to be an 
integer, so that the solutions are well-defined. 

Substitution of the expressions in (5.7) into (5.4)-(5.6) leads to the ordinary dif- 
ferential equation for R, 

R" + ^R' +(^k 2 - r ?pjR = 0, (5.8) 

where A; is a real constant, which is in essence the wavenumber. Satisfactrion of the 
surface boundary conditions leads to the dispersion relation 

u? = gktanhkh, (5.9) 



65 



which relates wave speed u/k to wavenumber, but is not of immediate concern here. 
The equation (5.8) can be written in the form 

(rR 1 )' + [k 2 r - R = 0, (5.10) 

and as such represents a Sturm- Liouville system of the form (5.1) (do not confuse the 
function r(x) in (5.1) with the polar radius r in (5.10)). 

The question of appropriate boundary conditions for (5.10) is an interesting one. 
We require 

R(a) = (5.11) 

to satisfy the no flow through condition on the tube wall, but the other condition 
required is more subtle. Since <fi satisfes Laplace's equation, we require V0 = 

^0 r ,-00,0 z ^ to be regular, and in particular, bounded. It follows from this that 
we require 

R(0) = 0, m^0, 

R'{0) = 0, m = 0. (5.12) 

The solutions of (5.8) are Bessel functions R = J m (kr), to which we will return in 
due course. 



5.1.2 Laplace's equation 



The second example we consider is that of solving Laplace's equation in spherical 
polars. This has a variety of applications, the most obvious of which are associated 
with properties of the Earth. For example, deviations of the Earth's gravity field 
from a purely radial one are usually expressed in terms of spherical harmonics, which 
are essentially the separable solutions of Laplace's equation in spherical coordinates. 
We suppose 

d 2 i> 



r 2 dr 



dr 



+ 



1 d 
r 2 sin 9 39 



sva.9 



(hp 

~d~e 



+ 



i 



r 2 sin 2 9 d(p 2 



0, 



(5.13) 



where (r, 9, </>) are spherical polar coordinates. Separable solutions exist in the form 

?/> = r n e im *e(9), (5.14) 



where we require m to be an integer in order that the solution be single- valued. Sub- 
stituting (5.14) into (5.13), we find that satisfies the ordinary differential equation 

-9 = 0. 



n(n + 1)6 + -A- [sin &]' - -% a 
sin 9 snr 9 



(5.15) 



It is convenient to write [i = cos 9, so that 

[(I-/' 2 )©'!' + 



n(n + 1) - 



1 d 



sin 9 d9 

2 



, and then 



m 



1-fi 2 



9 = 0. 



(5.16) 
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This is Legendre's equation. The solutions are usually denoted as P™(//), 
when m = 0, P° = P n are the Legendre polynomials, if n is an integer. 

We see that Legendre's equation is of Sturm-Liouville type, with A = n(n + 1). In 

spherical polars, V^> = (ip r: -ipe, -tps ), and thus the requirement that be 

\ r r sm / 

bounded implies that n = or n > 1, if the solution domain encloses the origin. If it 
encloses the poles at /J, = 1 and /j, = — 1, then we require 

6(±1) = if m = 0, 

9'(±1) = if (5.17) 



5.1.3 Instability in thermally activated shear flow 

A classical problem of some interest and relevance is that of thermal runaway. In 
its origins, this problem concerns combustion (why does a match light when you 
strike it?), but a different application occurs in the flow of fluids with temperature- 
dependent viscosity. The process of injection moulding (of, for example, telephones) 
fills a mould with a molten polymer, which then solidifies to form the manufactured 
object. If the flow rate of the polymer is high enough, charring can occur at the edges, 
which is evidently not desirable! 1 

The mechanism for this charring is that there is frictional heat released by the 
viscous flow of the melt, and this depends on the temperature, through the dependence 
of viscosity on temperature. The heat released increases the temperature, which 
further increases the heat released. This kind of circular reinforcement is called a 
positive feedback, and provides a mechanism for instability and, in some cases, what 
is called blow up: this is discussed further in chapter 9. 

A simple model which represents the injection moulding process is that of an 
axial flow of an incompressible fluid in a cylindrical tube. We use cylindrical polar 
coordinates (r, 9, z), and suppose that the tube wall is at r = a. If there is a purely 
axial fully developed velocity u(r,t), then the only stresses are a shear stress r and 
the pressure p, and these satisfy the axial component of the momentum equation, 
which is 

- £■ ^ 

Because of the purely radially dependent velocity, mass conservation (incompressibil- 
ity) is identically satisfied, and the acceleration terms in the momentum equation are 
identically zero. This latter situation is also approximately true for low Reynolds 
number (slow) flow. 

The viscous flow law is 

where fj, is the viscosity. We assume that \i = £t(T), where T is the temperature, 
whose evolution must therefore also be determined. The energy equation takes the 



J This was the subject of an industrial problem brought to Cornell University in 1975 by the Bell 
telephone company. 
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form (assuming axial symmetry) 



pc p 



OT OT 



k 



1 ( 0T\ 2 T 
r Or \ Or J Oz 2 



+ 



T 

p 



(5.20) 



where the terms represent successively heat advection, heat conduction and viscous 
dissipative heating: p is the fluid density, c p its specific heat, and k is the thermal 
conductivity. The last term in (5.20) is normally completely insignificant in ordinary 
laboratory scale fluid flows, but is important in the flow of molten polymers and 
plastics, and also in some larger geophysical flows. 

If we suppose that the temperature field is also fully developed, i. e., T = T(r, t), 
then the z derivative terms in (5.20) may be ignored. From the momentum equations 

0), and 



Op Op 
it follows that p' = — is constant (since also radial momentum implies — 
Oz Or 
therefore 

l / l / f lrdr 

r = ^rp', u = -2PJ r —> 

and the temperature satisfies the equation 



pc p 



OT _k ( OT' 
Ot r Or 1 Or 



+ 



r 2 p' 2 
4/i ' 



(5.21) 



(5.22) 



and is evidently uncoupled from the determination of the velocity field. We suppose 
temperature is prescribed on the tube wall, thus 



T = T on 



a. 



(5.23) 



Non-dimensionalisation 

We suppose that po is a typical value of the viscosity, and AT is the temperature over 
which it varies significantly; generally p decreases with T (think of maple syrup), and 
to be specific we will suppose 



p = Po exp 



(T - Tp) 
AT 



(5.24) 



We non-dimensionalise the temperature equation and boundary conditions by writing 

a 2 r 

r = a£, T = T + (AT)9, t = , (5.25) 



K 



where k = k/pc p is the thermal diffusivity. This leads to the dimensionless differential 
equation 



09 _ 1 

Ot~ Wt, 



.06 



2„0 



(5.26) 



with the boundary conditions 



on r = 1, 
at r = 0, 



(5.27) 
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this last condition following from the requirement that 9 be bounded at r — 0. The 
single dimensionless parameter a is given by 

a = , P (5.28) 

and we suppose that viscous heating is significant, i.e., a > 0(1). 

Steady state 

The steady state solution 0o(r) of (5.26) satisfies the nonlinear ordinary differential 
equation 

(W + +^ do = 0, (5.29) 



with boundary conditions 
We put 

so that ip satisfies 



O (1) = 0, 0' O (O) = 0. (5.30) 
= -ln£, 6 = 4z + i>, (5.31) 



iP" + ae^ = 0, (5.32) 
with boundary conditions 

ip = at z = 0, ip' — )• —4 as z oo. (5.33) 

A first integral of this can be found, which is 

\ip' 2 + ae^ = 8, (5.34) 

and the problem reduces to a quadrature, which determines the solution as 

[a 



if) = -2 In 
where c is given by 



cosh(22; + c) 

8 



(5.35) 



/ 8 

c = ±cosh _1 W-. (5.36) 
V a 

There are thus two solutions if a < 2 a/2, and none if a > 2y/2. For a greater than 
this critical value, the solution to the initial value problem for (5.22) blows up in finite 
time, and no steady state solution exists. It is this phenomenon which is associated 
with the charring which is observed if the injection moulding flow is too fast. 
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Linear stability 

We study the stability of a steady state solution by putting 

e = 9 + e. 



(5.37) 



Linear stability follows from assuming that (remember the variables are dimension- 
less) 9<1, and linearising the resulting evolution equation for O, to obtain 



e 7 



ld_ 



.50' 



+ a£ 2 e e °e. 



This has separable solutions of the form 

9 = e-P(0, 

where P satisfies 



together with boundary conditions 

p'(0) = 0, P(l) = 



P = 0, 



(5.38) 

(5.39) 
(5.40) 

(5.41) 



This is of the same form as the Sturm-Liouville problem (5.1). If written in terms of 
z given in (5.31), the problem becomes 



P" + [-ae~ 2z + 8sech 2 (2z + c)\ P = 0, 
where now P' = dP/dz; this is to be solved on < z < oo, with 

P(0) = 0, P'(oo) = 0. 



(5.42) 



(5.43) 



5.2 Sturm-Liouville theory 

The general Sturm-Liouville system is written in the form 

(py 1 )' + (Ar - q)y = 0, (5.44) 

where p, q and r are functions of the space variable x, and to be specific we solve the 
equation in the domain [a, b] and apply boundary conditions 

y(a) = y(b) = 0. (5.45) 

This is an eigenvalue problem: as we shall see, solutions exist only for a discrete 
set of eigenvalues Ai, A 2 , . . . , for which the corresponding functions </>i(x), <p2(x), . . . , 
are called eigenfunctions. We wish to demonstrate this, and elucidate some of the 
properties of these eigenpairs. 

We will assume that p, q and r are positive and smooth, that is to say continuously 
differentiable (e C 1 ), which provides the simplest case. It is possible to relax some of 
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these assumptions; in particular, it is often found (as in our examples (5.10), (5.16) 
and (5.40)) that p = at an end point; in this case the Sturm-Liouville system is 
referred to as a singular system. Depending on the singularity, the basic results which 
we describe below will still hold, with some modifications. 

The system (5.44) and (5.45) is a boundary value problem, which is why it is 
difficult to prove existence results. We therefore begin with a corresponding initial 
value problem, that is, (5.44) together with the initial conditions 



y(a) = 0, y'(a) = l. 



(5.46) 



We define z = py' and u 



where the matrix M is 



and 




so that 



u 



Mu, 



M 





\q-Xr 



P 




u(a) 




(5.47) 



(5.48) 



(5.49) 



By our assumptions, the matrix M is smooth and therefore Mu is Lipschitz contin- 
uous: hence Picard's theorem implies that there exists a solution which is unique so 
long as u is bounded. Being a linear system (with bounded coefficients), the solution 
can grow no faster than exponentially, and therefore the solution will exist for all x. 

The matrix M depends on A, so that the solutions do also, and in particular 
y = y(x, A). The dependence on A is smooth, and indeed analytic, as can be shown 

du 

simply by constructing the solution — of the differentiated system 



u' = Mu x + 





— rz 



(5.50) 



using the fundamental solution as in chapter 4. Thus the value of y at x = b, 

y(b,X) = y b (X) (5.51) 

is an analytic function of A. If y = at x = b (as well as x = a), i. e., y&(A) = 0, then 
we say A is an eigenvalue of the Sturm-Liouville system. The solution function y(x) 
is called the corresponding eigenf unction. 



5.2.1 Eigenvalues and eigenfunctions 

It will not have escaped the observant reader that Sturm-Liouville systems represent 
a generalisation of the simple harmonic oscillator 

y" + Xy = 0, y(0) = y(l) = 0, (5.52) 
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where p = r = 1, q = 0. Evidently the eigensolutions axe y n = sinn7rx, where n is an 
integer, and the corresponding eigenvalues are X n = n 2 ir 2 . We note that the eigenval- 
ues form a denumerable sequence which tends to infinity, and the eigenfunctions are 

orthogonal, in the sense that / y m y n dx = for m ^ n. It is as a consequence of this 

Jo 

that the functions {y n } form a basis for the expansion of functions in Fourier series, 
which is in fact complete, i. e., any smooth (or in fact piecewise continuous) function 
can be expanded in a Fourier series. We will find that these results carry across to 
the general Sturm-Liouville system. 

dy 

Let us denote w = — — . We then have 
o\ 

(py')' + (Xr - q)y = 0, 

(pw 1 )' + (Xr — q)w = —ry. (5.53) 
Multiplying by w and y respectively, subtracting and integrating, we find 

\p(y'w - yw')] b a = [" ry 2 dx > 0, (5.54) 

J a 

assuming r > 0. Since y = at a, b, and w = at a, it follows that 

{py'w)\ b = / ry 2 dx > 0, (5.55) 

J a 



dy 

and therefore w = — — / at x = b. It follows from this that the zeros of j/&(A) are 
isolated, since y^ is analytic, and if there an infinite number, the only limit point can 
be at ±oo. 



An eigenvalue exists 

It is by no means obvious that any eigenvalues will exist. To show that there is 
at least one, we consider the initial value problem (5.44) with the initial conditions 
(5.46). If p > 0, then we can suppose without loss of generality that p = 1, and we 
suppose that q and r are positive. 

If A < 0, then y" > as long as y > 0. Since y > for x close to a, this shows that 
y will continue to increase if A < 0, and therefore y&(A) > for A < 0. By continuity, 
this can be extended to state that y&(A) > for A positive and sufficiently small. 

Next we consider two solutions y\ and 2/2 with different values of A = Ai and 
A = A 2 , and we suppose Ai < A 2 . Multiplying the equations for y 1 and y 2 by y 2 and 
yi, respectively, and subtracting, we find 

(2/i2/2 - 2/12/2)' = ( A 2 - Xi)ryiy 2 . (5.56) 
Integrating, and then dividing by y\, we obtain 

(A2-A1)/ ry x y 2 dx 

= ^ > 0. (5.57) 

2/2 
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By l'Hopital's rule, we have — = 1 at x = a, therefore y\ > y 2 as long as y 2 > 0. 

2/2 

This implies that y&(A) is monotonically decreasing with A as long as it is positive, 
and also that, if xi(A) is the first zero of y (if it exists), then x\ is monotonically 
decreasing with A; in particular, x\ > b for small positive A. 

On the other hand, it is fairly obvious that if A ^> 1, then by comparison of the 
solution of (5.44) (with p = 1) with that of 

z" + fiz = 0, z(a) = 0, z'(a) = 1, (5.58) 
where A > ^ max ^ an d ^ [ s positive, we can show that a?i(A) does indeed exist, and 

^min 

xi — >• as A — > oo. Since x\ varies continuously with A, is monotonically decreasing 
with A, and x\ > b for sufficiently small A, this implies that there is a unique positive 
value Ai where xi(Ai) = b, and this is an eigenvalue of the Sturm-Liouville system. 



Extension to the case p ^ 1 is trivial, by defining the new space variable X = 




and suitable redefinition of q and r. 

This proves that there is an eigenvalue. Moreover, we have in passing proved that 
all the eigenvalues are positive, their only limit point is oo, and that the minimum 
eigenvalue Ai has a corresponding eigenfunction 4>i(x), say, which is of one sign. 
In addition, it is simple to repeat the argument for A > Ai (simply redefine the 
left endpoint a\ to be at x\{X)) to show that there must be an infinite sequence of 
eigenvalues, and the sequence is denumerable since it can be exhaustively constructed 
in the above way. 



Eigenfunctions are orthogonal 

A fundamentally important property of the eigenfunctions of Sturm-Liouville systems 
is that they are orthogonal. Orthogonality is a property which we associate most easily 
with vectors. In three dimensions, two vectors a and b are at right angles if their scalar 
product ^2 ciibi is zero, and we extend the definition to less geometrically accessible 

i 

n-dimensional vectors simply by allowing the same definition of scalar product. 

The language we are using to discuss eigenvalues is very suggestive of the fact 
that we are here dealing with the same algebraic and geometric structures as occur 
in finite-dimensional vector space. There is, in fact, a direct analogy: functions are 
essentially (countably) infinite dimensional vectors; the finite number of eigenvalues 
of a matrix is replaced by the denumerable sequence we have discussed, and just as 
the eigenvalues of a real symmetric (or self-adjoint) matrix are real and orthogonal, 
so also a Sturm-Liouville system is self-adjoint (this will be elaborated in chapter 6), 
and has real eigenvalues which are orthogonal. The definition of a scalar product 
for a function space is called an inner product, and its definition is the extension of 
the finite sum ^ a;6; to its limit — an integral. To be specific, we define the inner 

i 

product of two functions u and v with (positive) weight function r(x) to be 

rb 

(u,v) = / ruvdx, (5.59) 
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and we say that u and v are orthogonal if (u,v) = 0. 

For the eigenf unctions y 1 and y 2 with corresponding eigenvalues Ai and A 2 , we 
find from the differential equation that 

WiV2 ~ 2/12/2)]' + (^1 - M)rym = 0. (5.60) 
Integrating from a to b and using the boundary conditions, we find 

rb 

/ ryiy 2 dx = 0, (5.61) 

J a 

that is, yi and j/2 are orthogonal. 

Just as we construct bases for vector spaces from sequences of orthogonal vectors, 
so we can attempt to construct 'bases' of orthogonal functions for function spaces. 
This is the realm of Fourier series. We have an infinite sequence of eigenfunctions, and 
if we take these to be orthonormal, in the sense that (</>„,</>„) = 1, then a candidate 
Fourier expansion of an arbitrary function f(x) is 

/ = (5-62) 

i 

and the issue of whether this candidate expansion is actually correct resides in whether 
the eigenfunctions form a 'complete' sequence. The answer is that they do, but a 
complete proof of this lies somewhat beyond the scope of these notes. We will sketch 
a proof in chapter 6. 



A variational principle 

Multiplying the Sturm- Liouville equation for y by y and integrating, we find that 
each eigenvalue A satisfies 

rb 

/ {qy' 2 +py 2 )dx 

A = — T£ • (5.63) 

/ ry 2 dx 

J a 

Since this expression is homogeneous, we can take the denominator to be equal to 
one without loss of generality. In different phrasing, we define the norm of a function 
u to be 

1/2 



\u\ 



!b ~) ' 

f ru 2 dx \ , (5.64) 



that is, \\u\\ = y(m,m), in analogy with the Euclidean vector norm, so that we may 
take 1 1 y \ \ = 1 in (5.63). In that case 

A= [ b (qy ,2 + P y 2 )dx. (5.65) 

J a 

It is a simple observation of the calculus of variations that the functional X[y] de- 
fined by (5.65) (or (5.63)) takes stationary values when the Sturm-Liouville equation 
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is satisfied, and in fact is also a local minimum (since p and q are positive). The first 
eigenvalue Ai is equal to the global minimum of X[y], and then the next one A 2 is the 
global minimum of X[y] subject to (y,yi) = 0, and so on. It is somewhat easier to 
prove these statements when the Sturm-Liouville system is rewritten as an integral 
equation, and this is done in chapter 6. 

5.2.2 Singular Sturm-Liouville systems 
5.3 Comparison methods 

Related to eigenvalue problems are the oscillation and comparison theorems, which 
apply generally to second order differential equations, and concern the existence and 
location of zeros of the solutions. We will prove four related results of this type. 
Suppose first that 

(py[)'-Qiyi = 0, 

(py'2)' -Q2V2 = 0, (5.66) 

and q 1 > q 2 . Then y 2 has at least one zero between two successive zeros of y\. For, 
if this is not the case, let us take x\ and x 2 to be two successive zeros of yi, and 
suppose without loss of generality that y\ > for x\ < x < x 2 . If J/2 is non-zero in 
the same interval, we may suppose that it is also positive. Using a by now familiar 
construction, we have that 

\p(y'm - 1/11/2)]' = (ft - Qi)ym > o, (5.67) 

and integrating this between x\ and x 2 yields 

(my'i)\ X2 > (py2y[)\ xi , (5.68) 

which is impossible since we must have y[\ < and y[\ > 0. 
Next we show that the zeros of independent solutions of 

(py')' -qy = o (5.69) 

interlace. Suppose y\ and y 2 are independent solutions of (5.69): then their zeros are 
distinct (since if y\ = y 2 = at some value of x, then Picard's theorem implies that 
one is a multiple of the other). If x 1 and x 2 are successive zeros of y 1 , then as above, 
we can find that 

W1V2 - viy'X = 0, (5.70) 

whence 

W)L = (pwi/DU- ( 5 - 71 ) 

Because y[(xi) and y[(x 2 ) have opposite signs, it follows that so also do y 2 (xi) and 
y 2 (x 2 ), and thus that y 2 has a zero between x 1 and x 2 . The same argument works 
the other way round. This proves that the zeros of y 1 and y 2 interlace. 
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Now we turn to the Sturm- Liouville eigenvalue problem. If y\ and y 2 are two 
eigenfunctions with eigenvalues Ai and A 2 , so that 

(py[)' + (Air - q)y 1 = 0, 

(Py'2)' + (V - q)V2 = 0, (5.72) 

let us suppose that x\ and x 2 are two successive zeros of yi, with y\ positive between 
them, and suppose also that y 2 > in this interval. If Ai < A 2 , then 

\VV'MI\ = (*2 - Ai) H ry lV2 dx > 0, (5.73) 

J XI 

which is impossible. It follows that if Ai < A 2 , then the eigenfunction y 2 has a zero 
between any two zeros of y\ . 

Finally we prove again the result on the number of eigenvalues. For (5.44), we 
define 

*=r.m (5 - 74) 

so that with y' = dy/d£, we have 

y" + [Apr -pq]y = 0. (5.75) 
We compare the solutions of this equation with those of 

z" + (Ai? - Q)z = 0, (5.76) 
where R and Q are positive constants such that pr > R and pq < Q, and we suppose 

A > ^ > 0. Then Q — XR > pq — Apr, and the first result in this section (after (5.66)) 
R 

then implies that y has a zero between any two of z. Since by construction, z has an 
infinite number of zeros, it follows that so also does y. 

Denote these as x n (X), n = 1,2, . . .. From (5.55) it follows that y'y\ > at x n , 
where also y[x n (X), A] = 0. Thus we find 

<( A ) = -ff'< ' ( 5 - 77 ) 

and thus all the zeros decrease as A increases. From all this we can (again) deduce 
that there are an infinite number of eigenvalues X n , with X n — > oo as n — > oo. 



5.4 Notes and references 

A good source for the material on separation of variables and the derivation of the 
main equations of mathematical physics is the classic text by Jeffreys and Jeffreys 
(1946). This book has been reprinted many times and is still in print. 

There are many books on differential equations which cover Sturm- Liouville equa- 
tions. A classical text is that by Courant and Hilbert (1937). Other books, no less 
valuable, are those by Mackie (1965), Coddington and Levinson (1972), and Stakgold 
(2000). Oscillation and comparison theorems are nicely discussed in the little book 
by Burkill (1956). Simmons (1972) is a very accessible entry level text. 
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Exercises 



5.1 Show that the second order differential equation 



w" + a(x)w' + [b(x) + Xc(x)w] = 



can be written in the Sturm- Liouville form 



(py')' + (Xr - q)y = 0, 



by suitable definitions of the functions y, p, q and r. 
If the boundary conditions are 



w' + aw = at x = 0, 1, 



what are the corresponding boundary conditions for yl 



5.2 Bessel's equation of order v is 



w" + -w 1 + 1 - 




w = 0, 



and the solutions are called Bessel functions of order v. Show that the solutions 
of the Sturm-Liouville equation 



can be written as Bessel functions of order m. 

By seeking solutions near z = of the form w ~ z^, show that if v ^ 0, then 
independent solutions w± are possible with w± ~ z ±v . Deduce that if we require 
R(0) = when m is a positive integer in (*), then we must have R oc w + (kr). 

(w + is normally written J u , and is the Bessel function of the first kind; the 
Bessel function of the second kind is singular at z = and is denoted Y v .) 

Show that if v = 0, then independent solutions are possible in which 



5.3 A function u(x,t) satisifies the nonlinear reaction-diffusion equation 




R = (*) 




Y - 



\nz. 



Ut = u xx + f(u) on [0,1], 



with boundary conditions 



u = on x = 0, 1. 
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Suppose that u (x) is a steady solution. By writing u = u + U, derive a 
linearised equation for small perturbations U. Show that is has separable solu- 
tions of the form U = v(x)e at , and deduce that v satisfies the Sturm- Liouville 
equation 

v" + [s(x) - a]v = 0, v(0) = v(l) = 0, 

where 

s(x) = f'[uo(x)]. 

By consideration of a suitable variational principle, or otherwise, show that if 
s(x) < 0, then all eigenvalues a are negative, and hence the steady state is 
stable. 

If s(x) > 0, show, by writing q = s max — s and A = s max — a, that instability 
can occur (a > 0) if 

f (v 12 - s{x)v 2 )dx 

JO 

JO 



min ^ j < 0. 

° dx 



(This shows that the constraint that q > in Sturm-Liouville theory is not 
necessary.) 

5.4 A shear thinning fluid (like tomato ketchup or Guinness) has a viscosity which 
decreases with increasing stress. Suppose such a fluid flows in a cylindrical tube 
of radius a, and has viscosity 



(T-Tp) 
AT 



where r is the shear stress, and the power law exponent n > 1 (n = 1 corre- 
sponds to Newtonian flow). 

Show that a suitable dimensionless formulation of the temperature profile across 
the tube satisfies the dimensionless equation 



06 _ 1 d 



09 



with the boundary conditions 

= on r = 1, 

and 

Q r = at r = 0. 
Give the definition of the parameter j3. 

Find a steady solution of this problem, and find a Sturm-Liouville eigenvalue 
problem, whose eigenvalues determine the linear stability of the steady solution. 
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Chapter 6 



Integral equations and 
eigenfunctions 



6.1 Integral equations 



Like differential equations, the formulation and solution of integral equations is of 
interest in its own right. In this chapter, we will however focus on one particular 
class of integral equations, those which arise through the reformulation of a Sturm- 
Liouville differential equation system as a Fredholm integral equation. More generally, 
we begin by illustrating how this reformulation occurs. 

6.1.1 Volterra integral equations 

Volterra integral equations are those in which the integral occurs in indefinite form 
(i. e., the space variable x occurs as the upper limit of the integral). Typically Volterra 
integral equations are associated with initial value problems. Consider for example 
the inhomogeneous differential equation 



We can solve this using the method of variation of parameters, as follows. Suppose 
that D\ and y 2 are independent solutions of the homogeneous equation (i.e., with 
f(x) = 0) satisfying 



(py 1 )' -qy = - f, 



(6.1) 



together with the initial conditions 



y(a) = y'(a) = 0. 



(6.2) 



y 1 (a) = 0, y[(a) = l, 
1/2 (a) = 1, y' 2 (a) = 0. 



(6.3) 



We seek a solution in the form 



y = A(x)y 1 +B(x)y 2 , 



(6.4) 
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where the coefficient functions A and B are to be found. We require two equations 
for A and B, and we take the first of them to be 

A' yi +B'y 2 =0. (6.5) 

With this choice, we have that 

y' = Ay' 1 + By' 2 , (6.6) 

and the initial conditions for y in (6.2) require A and B to satisfy 

A(a) = B(a) = 0. (6.7) 

Multiplying (6.6) by p, differentiating, and substituting into (6.1), we find the 
second equation for A and B: 

y[A' + y' 2 B' = -L (6.8) 
Inverting the matrix equation formed from (6.5) and (6.8), we find 

(t)-Mi? )(-/)• 

where W = y\y' 2 — y[y2 is the Wronskian. Directly from the equations for y\ and y 2 , 
we know that pW = C is constant (non-zero, since the solutions are independent), 
and therefore the solutions for A and B satisfying the initial conditions (6.7) are 

A = ^[ a v2(zmz)dt, b = -± r yi (om)dt (6.10) 

We can thus write the solution compactly as 

y= [ X K(x,0f(0d{, (6.11) 

J a 

where the integral kernel is defined by 

K(x, = VMV*® -ViiOtoW . ( 6 . 12 ) 

This simply defines the solution. An integral equation occurs if we put / = Xry, 
thus regaining the Sturm- Liouville equation (with initial conditions) of chapter 5. In 
this case, (6.11) becomes 

y = X f K(x,0r(0y(0dt (6.13) 

J a 

This is now an integral equation of Volterra type, because of the upper limit x in the 
integral. In principle it can be solved using Picard iteration. 

More generally, if we seek eigenfunctions satisfying y(a) = 0, y(b) = 0, they can 
be formally constructed by solving the initial value problem with y(a) = 0, y'(a) = 1, 
thus 

y = X K(x,0r(0y(0dC + yi(x), (6.14) 

J a 

and then choosing A so that y(b) = 0. Practically, though, this is not a very good 
idea. 
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6.1.2 Fredholm integral equations 

Now we consider the boundary value problem for (6.1), i. e., 

Ly = ( P y')'-qy = -f, y(a) = y(b) = 0. (6.15) 

The solution for this can be written in terms of a Green's function, which is con- 
structed from independent solutions y\ and y 2 , which we now take to be defined as 
solutions of the homogeneous version of (6.15) (i. e., / = 0) with 

yi (a) = 0, y[(a) = l, 

y 2 (6) = 0, y' 2 (b) = l. (6.16) 

We explicitly assume that zero is not an eigenvalue for the corresponding Sturm- 
Liouville system, so that this choice gives independent solutions. 
The Green's function G(x,£) satisfies the homogeneous equation 

LG = (pGj -qG = in x < £, (6.17) 

with 

[Gt- = 0, \pG']l + _ = -l, (6.18) 

where the primes denote differentiation with respect to x, and indicates the 
jump from x just less than £ to x just greater than £. In terms of the independent 
homogeneous solutions y 1 and y 2 , the Green's function is explicitly given by 

g <*-«=< «(«&(.)' (6 - i9 > 

J X > S) 



c 

where C = pW is constant (and non-zero), W being the Wronskian for this pair yi, 

If now we consider a Sturm- Liouville system in which / = \ry, then the eigen- 
functions y and eigenvalues A satisfy the Fredholm integral equation 

y(x) = X [ b G(x,0r(0y(0dt (6.20) 

J a 

This is called a Fredholm integral equation because the limits on the integral are 
fixed. 

We may also derive an inhomogeneous Fredholm equation by consideration of the 
inhomogeneous Sturm-Liouville equation 

(py>y + (\r-q)y = -f, (6.21) 

subject to y(a) = y(b) = 0. This leads to the integral equation 

y(x) = \[ b G(x,Z)r(Z)y(Z)d£ + g(x), (6.22) 

J a 



where 



g(x)= [ b G(x,0f(0dt (6.23) 

J a 
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6.1.3 The delta function 



Although not essential to the discussion, it is difficult to go past the derivation of 
the Green's function without some reference to the delta function 5(x). The delta 
function is like an administrative tool; it makes things easier to treat, once you have 
the confidence to use it, but at first sight it is very strange. 

For a start it is not a function. It is called a generalised function. It has the 
property that it is zero everywhere except at x = 0, where it is 'infinite enough' that 

/oo 
S(x)dx = l. (6.24) 
-oo 

Obviously this is not really satisfactory. There are two formal ways of thinking 
about generalised functions. The first is that we think of a generalised function as 
being defined as a 'limit' (more precisely, an equivalence class) of a family of smooth 
functions; for example, the delta function is associated with the family of functions 

e~ nx as n — > oo. (Note that there are many such families.) 

But the main way of manipulating generalised functions is not to manipulate them 
directly at all, but to define their action by means of integrals. If {g n (x)} is a family 
of functions in the equivalence class of a generalised function g(x), then we define the 
integral 

/OO POO 
f(x)g(x) dx = lim / f(x)g n (x)dx. (6.25) 
-oo n ~ >OC J — OO 

Using this definition, it is possible to prove that one can manipulate generalised 
functions directly by using this integral property, so that for example, one can define 
the delta function through the requirement that 

/oo 
f(x)5(x)dx = f(0) (6.26) 
-oo 

for all smooth functions f(x). 

With this rudimentary description of the delta function in hand, we can more 
conveniently define the Green's function G to be the solution of 

(PC')' -qG= -5{x - (6.27) 

which satisfies the boundary conditions G = at x = a and x = b. Integration of 
(6.27) across x = £ then implies the second equality in (6.18). 

6.1.4 The Fredholm alternative 

The Fredholm alternative concerns the solution of the inhomogenous Fredholm equa- 
tion (6.22). To put this in more familiar shape, it is convenient to define an integral 
operator 

Qu= [ b G(x,0r(0<0dt (6.28) 

J a 

We should think of this in the same way that we (should) do of functions and linear 
transformations. Both define mappings: a (real- valued) function of one (real) variable 
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f(x) defines a mapping from the real numbers to the real numbers; a linear transfor- 
mation M on a finite-dimensional vector space maps a vector v to another vector Mv. 
In a similar way, an integral operator such as Q in (6.28) or a differential operator 
such as L in (6.15) maps a function to another function. These functions reside in a 
suitable function space, and the issue of what a suitable function space should be is 
of technical concern, and forms a focal point of functional analysis. It will only be 
of diversionary interest in these notes, but it is as well to be aware that many of the 
proofs of the subject really require attention to the detail of what function space the 
operators are defined on. Unlike, say, linear transformations, it is usual for mappings 
of functions to be defined on one function space, but to have their range, or image, 
in another function space. As a simple example, we might define the operator L in 
(6.15) as acting on the space C 2 [a, b] of twice continuously differentiate functions on 
[a, b], but its image is then clearly only the continuous functions C[a, b]. Apart from 
this, we should think of operators on function spaces in much the same way as we do 
linear (or nonlinear) transformations on vector spaces. 

With the definition of Q in (6.28), we might take u € C[a, b], and then its image 
is also in C[a, b] if, as we suppose, G is continuous. The inhomogeneous Fredholm 
integral equation (6.22) can then be written in the compact form 



The Fredholm alternative is essentially the same statement we would make about 
this equation if y and g were finite-dimensional vectors, and Q were a matrix. It 
is this: either the homogeneous equation g = has no non-zero solution, and then 
(6.29) has a unique solution; or, the homogeneous solution has a solution <j>(x) (an 
eigenfunction), and then the inhomogeneous equation has no solution at all unless 
g is orthogonal to (p, in which case there are an infinite number of solutions (a one 
parameter family). The particular form of the statement above applies to self-adjoint 
linear operators Q\ this was touched on briefly, as was the definition of orthogonality, 
in chapter 5; it is discussed further in section 6.2.2 below. The proof of the Fredholm 
alternative will emerge as we discuss eigenfunctions further. 

6.2 Eigenfunction expansions 
6.2.1 Neumann series 

Just as Picard's iterative method for Volterra integral equations leads to a convergent 
series form of solution, so also does an iterative method for the inhomogeneous Fred- 
holm integral equation. It is most easily written using the operator notation, thus we 
solve (6.29) by defining the sequence of functions 



y = XGy + g. 



(6.29) 



VO =9, Vn = XQVn-l + g- 



(6.30) 
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This leads to the sequence y n = ^ X n Q n g, and suggests that we may write the solution 

o 

as the formal series 

oo 

y = J2\ n G n g. (6.31) 
o 

This is known as a Neumann series, or Neumann expansion. Note that the solution 
can be written in the formal way 

y=(I-\g)- 1 g, (6.32) 

where / is the identity operator, and that the Neumann series results from the formal 
expansion of (/ — XQ) -1 as a Taylor series. In fact, the issue of convergence of this 
series (and hence of validity of this formal solution) requires consideration of the norm 
of the operator Q. So now we revisit the definitions given in chapter 5 of norms, and 
we extend their definition from norms of functions to norms of operators. 



6.2.2 Inner products and norms 

We recall the definition of an inner product from chapter 5. If r(x) > 0, then we can 
define an inner product on a suitable space of functions as 

(u,v) = / ruvdx. (6.33) 

J a 

In reality, the space is often taken to be L 2 , the space of Lebesgue square-integrable 
functions, but in practice, we will limit ourselves to considering the space of contin- 
uous functions. In general, an inner product (u,v) is a quantity which is symmet- 
ric ((u,v) = (v,u)), bilinear (i.e., linear in both u and v arguments), and satisfies 
(u,u) > 0, with equality only if u = 0. 1 

We have already mentioned the idea that an operator be self-adjoint. The adjoint 
L* of an operator L is defined (uniquely) by the property that 

(Lu,v) = (u,L*v), (6.34) 

and the operator is self-adjoint if L = L*. It is easy to see that the Sturm- Liouville op- 
erators of chapter 5 and the integral operators of the present chapter are self-adjoint. 
It is this property which ensures the eigenvalues are real, and the eigenfunctions are 
orthogonal. 

A space which has an inner product is called an inner product space, and an inner 
product space is also endowed with a norm: 

| |u|| = {u,uf /2 . (6.35) 

In general, a norm \ \u\\ is a distance- like measure of an element u of a space which is 
positive (unless u = 0), linear under scalar multiplication, and satisfies the triangle 
inequality 

||w + < |H| + |H| V (6.36) 

1 As we do throughout this discussion, we ignore the possibility that u and v are complex, and 
restrict ourselves to the case where they are real. 
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Spaces having norms are called metric spaces, and they allow the notion of conver- 
gence. If every convergent sequence of elements converges to an element of the same 
space, then the space is called complete. A complete, normed linear space B (i.e., if 
u e B and v e B, then au + (3v e B for all ct,/3 € R) is called a Banach space. A 
complete inner product space is called a Hilbert space. 

In an inner product space, we have that ||m + sv\\ 2 = \\u\\ 2 + 2s(u,v) + s 2 |H| 2 is 
always positive, and from this follows the Cauchy-Schwarz inequality: 

(u,v) < \\u\\ \\v\\. (6.37) 

It is as a consequence of this that the inner product defines an associated norm (6.35), 
since the triangle inequality can be deduced from the Cauchy-Schwarz inequality. 

One of the reasons we bother so much with what space a function is in is that 
we want convergent sequences to converge in the same space: we want to deal with 
complete spaces. It is not difficult to see that apparently nice looking spaces, such as 
the continuous functions on [0, 1], denoted C[0, 1], are not complete. Endowed with 

Ui 1 1/2 

u 2 dx> , the sequence of functions u n = e~ nx converges, in the 

sense that ||it„|| — > as n — > oo. But the limit function is clearly u(0) = 1, u(x) = 0, 
x > 0, which is not continuous. So C[a, b] is not complete with this norm. The reason 
that the square-integrable function space L 2 is of interest is that it is complete with 
this norm. 

Given the property of a norm, we can also define norms of operators. We define 
the norm of an operator Q to be 

||0|| = sup ^f, (6.38) 

IMI 

or equivalently, \\Q\\ = sup ||(?it||. We leave it as an exercise to show that this 

Nl=i 

definition indeed defines a norm. It is then immediate that the Neumann series 
(6.31) converges if 

|A| < jjlf, (6.39) 

and is an analytic function (or functional) of A within this circle of convergence. 
The series convergence is compromised by a singularity on the circle of convergence, 
and this singularity occurs at an eigenvalue. More generally, analytic continuation of 
(I — At/) -1 to other complex A is limited only by singularities at the (real) eigenvalues 
Ai, A 2 , . . .. 

It remains to show that the integral operator Q defined by (6.28) is indeed bounded, 
so that its norm exists. For each x, we have in fact 

Qu = (G,u) < \\G\\ \\u\\, (6.40) 

and thus 

< [ b [ b r(Or(v)G 2 (^v)^drj, (6.41) 

J a J a 



\\Qu 


1 1 2 




\u\\ 


2 
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whence Q is bounded if the right hand side of (6.41) is bounded (for example, if G is 
continuous), and then 

J a r(Or(v)G 2 (tv)d^d V \ . (6.42) 
This property defines Q as a Hilbert- Schmidt operator. 

6.2.3 Spectral theory 

The discussion of eigenvalues and eigenfunctions sits within the more general func- 
tional analytic subject of spectral theory. For a general operator Q, the values of A 
where the equation (/ — XQ)u = has a non-trivial solution are the eigenvalues, and 
in the present situation, where we have a Hilbert- Schmidt operator, the eigenvalues 
constitute the entire spectrum. More generally, the eigenvalues constitute the point 
spectrum. If we suppose that Q maps functions in a space U to U (e. g., we might have 
U = C[a, &]), then values of A for which (J — XQ)u = has no non-trivial solution, 
but where the range of (I — XQ) is not the whole of U, constitute the continuous 
and residual spectrum; the continuous spectrum is that part for which (J — XQ) -1 is 
unbounded. 



6.3 Variational methods 

Most of what we have said so far applies to general Fredholm integral equations with 
real symmetric kernels, and this will continue to be the case. For the particular 
integral operators which are the inverse of the Sturm-Liouville operators, we already 
know that the eigenvalues are characterised by a variational principle, specifically 

Ai = min y\qu 2 + pu' 2 ) d£ J , (6.43) 

with similar variational principles for the other eigenvalues. It is easiest to prove such 
assertions when the problem is formulated as an integral equation. 

6.3.1 A variational principle 

Now let us address the Fredholm integral eigenvalue problem 

rb qj 

Gu = j a G(x,Or(CHOdC=y (6.44) 

where we assume G is real and symmetric, so that Q is self-adjoint. The Cauchy- 

Schwarz inequality states that (Qu, u) < \ \Qu\ \ \ \u\ \ < \ \Q\ \ if | |it| | = 1, and it is simple 

u 

to show that when this upper bound is obtained, then Qu = — , and thus we have 

Ai 

sup (gu,u) = ^ = \\g\\. (6.45) 

IHI=i Ai 
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All that is missing from this characterisation of the first eigenvalue is a proof that 
the supremum is actually attained. This takes us to the periphery of our exposition. 
Prom the definition of a supremum, we have that there is a sequence of functions 
{u m } such that (Qu m ,u m ) — > \\Q\\ = 1/Ai as m — > oo. For u m € U = C[a,b], Qu m 
is an equicontinuous, uniformly bounded sequence of functions (an equicontinuous 
family is one which is uniformly continuous over both x and m). The Arzela-Ascoli 
theorem states that an infinite equicontinuous, uniformly bounded sequence of func- 
tions contains a convergent subsequence. Taking this sequence to be {u m }, this shows 
that Qu m — > (pi(x), where </>i is continuous, and it is straightforward to show that 

G<pi = — -; 0i is the first eigenfunction. 
Ai 

Now consider the modified integral equation 

rb qj 

g 2 u = J a G 2 (x,Or(0«(0« = y (6-46) 

where 

G 2 (x,a = G(x, 5 )-*MM). (6.47) 
As in the process of Gram-Schmidt orthogonalisation, we can write 

u = v + a<J> 1 , (6.48) 

where 

v = u - (tt,0i)0i, a=(u,(j) 1 ), (6.49) 

and thus (f,0i) = 0. By direct calculation, we find that (G2U,u) = (Qv,v), and as 
before, we can deduce that 

sup (Q2U, u) = sup (Qv, v) = (6.50) 

IHI=1 («,«i)=0 A 2 

is obtained at a second eigenfunction 2 , which is orthogonal to the first. 

We can go on iteratively in this way to construct a sequence of operator equations 

rb 111 

g n u = J a G n (x,0r(0<0d£ = (6.51) 

where 

G n (x,0 = G(x,0 - £ (6.52) 

1 A? 

and we find that the n-th eigenvalue X n satisfies 

^- = sup (G n u,u) = \\Gn\\, (6-53) 

||«||=1 

and this bound is attained when u = (j) n , the n-th eigenfunction. The eigenfunctions 
(with \ \4>j\ \ = 1) form a denumerable orthonormal set. 
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We already know from chapter 5 that X n — > oo as n — > oo, at least for the Green's 
function kernel G(x,£). We can also show this directly by calculating 

ff G 2 n (x, 0r(s)r(O dx d£ = ff G 2 (x, 0r(s)r(O d£ -"f^ (6.54) 

Ja ./a J a J a ^ 

from which it follows that the infinite sum converges, and we have a form of Bessel's 
inequality: 

Ei2 < / / G 2 (x,Or(x)r(OdxdC. (6.55) 

l Aj JaJa 

In particular, there is no finite limit point for Aj, and Aj — > oo as j — > oo. 

Since ||(7„|| = 1/A„ — > as n — > oo, it follows that for any m on the unit sphere 
= 1, the sequence — > as n — >• oo. As before, Arzela-Ascoli implies 
that there is a subsequence such that {G n u} —> 0, and since Q n is equicontinuous, 
this implies that Q n — > 0, and thus that G n — >• 0; hence we have the expansion of the 
kernel in the form 



oo 



G(x,0 = J2 YiiX l m - (6-56) 
i A? 



6.3.2 Fourier series 



Having constructed an orthonormal sequence of eigenf unctions, we would like to be 
able to solve some of the general initial value problems which give rise to the Sturm- 
Li ouville problems of chapter 5. Whether we can do this depends on whether the 
orthonormal sequence of eigenfunctions is complete. To see this, let us suppose that 
the function f(x) € C 2 [a, b] (this is not actually necessary, but makes the demonstra- 
tion simpler). With the Sturm-Liouville differential operator defined by (6.15), we 
have that u = —Lf is continuous, and its inverse is / = Qu. We wish to show that / 
can be expanded in a Fourier series, 

oo 

/ = £(/,fc)fc- ( 6 - 57 ) 

1 

Note that (/, <pj) = (Qu, fy) = (u, Q4>j) by the property of self-adjointness, and this 
equals —(u,(j)j). Hence we can deduce that 



A 3 



m—l m—1 



f ~ £ (/, 4>3)<f>3 = $\- U - E ( U ' faWj] = GrnU, (6.58) 
1 1 

and thus 

m ~ 1 I lull 

II/- £(/,&)&! I = \\<3™v\\ < \\0m\\\\u\\ = ^ ^0 as m^oo. (6.59) 

A r , 



1 



88 



Again, equicontinuity of {Gm} implies that / — J2T 1 (f, (j)j)(f)j converges to zero, and 
hence (since / is continuous) we have 

oo 

/ = £(/>fc)fc- (6-60) 
1 

As a corollary (although it is in fact equivalent), we have Parseval's relation: 

oo 

ll/H 2 = £l(/^)| 2 - (6-61) 
l 

A version of (6.61) with inequality (left hand side > right hand side) is Bessel's 
inequality. 

6.3.3 The Fredholm alternative revisited 

6.3.4 Rayleigh-Ritz method 

6.4 Notes and references 

There are a good number of books on integral equations. Most of the material in 
this chapter can be found in Courant and Hilbert (1937), a classic exposition, or 
Stakgold (2000), which has a more up to date flavour. Generalised functions are 
treated elegantly in the little book by Lighthill (1958). 

Exercises 

6.1 Show that the Wronskian W = yiy' 2 — y[y2 of two solutions yi, y<i of the differ- 
ential equation 

\p(x)y']' - q(x)y = 0, p(x) > 0, 

satisfies pW = C, a constant. Show that the solutions are dependent (one a 
multiple of the other) if C = 0, and independent if C ^ 0. 

Use the method of variation of parameters to find the solution in integral form 
of 

y" + cu 2 y = -f(x), y(0) = y'(0) = 0. 
By constructing the Green's function, find the solution of 
y" + u 2 y=-f(x), y (0)=y(l) = 0, 

in the form 

y = /"'coco/code- 

Jo 
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6.2 The integral operator Q is defined on the space of continuous functions on [a, b], 
C[a,b], by 

Qu= /*G(z,Or(0«(flde. 

J a 

where G is real, symmetric and continuous on [a, b] x [a, b], and r is positive and 
continuous on [a, b]. 

Suppose that the eigenfunctions <p n with corresponding eigenvalues X n of 

(/ _ \g)<p = o 

form a complete, orthonormal sequence, and that a continuous function / has 
a Fourier expansion 

i 

By writing (/ — At/) -1 = ^ X n Q n , show that a formal solution of the inhomo- 

n 

geneous Fredholm integral equation 

(/ - XG)y = f 

can be found in the form 

EQ A j (pi . . 

Show that the derivation is only valid for |A| < jj^yp ^ na ^ ^ ne ser i es m (*) 
converges for all A except the eigenvalues (where there are simple poles). 

6.3 In deriving the Neumann series 

oo 


for the solution of 

(/ - \Q)y = g, 

we implicitly assumed that the operator Q is associative, that is, QoQ n = Q n oQ 
for all positive integral n. Demonstrate this explicitly, first for n = 2, and then 
by induction for all n, by deriving an expression for the kernel of the integral 
operator 

J a 

6.4 An integral operator is defined by 

Gu= l*G{z,t)T(tMt)dZ, 
90 



where G is real, symmetric and continuous, r is positive and continuous, and u 
is continuous. 



Show that 



rb 



1/2 



= i I ru 2 d£ 



defines a norm for u, and that 

\\Q\\ = sup 

defines a norm for Q. Hence show that 

rb r b 



\Qu\ 
\\u\\ 




1/2 



\\G\\< / / r(Or( V )G 2 ^, V )d^drj 



If G is the Green's function calculated for Bessel's equation of order zero on the 
interval [0, 1], is bounded? (See also question 5.2.) 

By consideration of the kernel 



i A i 



or otherwise, derive the Bessel inequality 



E^2 < / / G 2 (x,0r(x)r(0dxdt, 



a J a 



and deduce that \ n — > 00 as n — > 00. 
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Chapter 7 
Waves and shocks 



Any introductory course on partial differential equations will provide the classifica- 
tion of second order partial differential equations into the three categories: elliptic, 
parabolic, hyperbolic; and one also finds the three simple representatives of these: 
Laplace's equation V 2 w = 0, governing steady state temperature distribution (for 
example); the heat equation u t = V 2 m, which describes diffusion of heat (or solute); 
and the wave equation u u = V 2 «, which describes the oscillations of a string or of a 
drum. These equations are of fundamental importance, as they describe diffusion or 
wave propagation in many other physical processes, but they are also linear equations; 
however the way in which they behave carries across to nonlinear equations, but of 
course nonlinear equations have other behaviours as well. 

In the linear wave equation (in one dimension, describing waves on strings) u u = 
c 2 u xx , the general solution is u = f(x+ct)+g(x — ct), and represents the superposition 
of two travelling waves of speed c. In more than one space dimension, the equivalent 
model is u tt = c 2 V 2 m, and the solutions are functions of (k.x ± ut), where u is 
frequency and k is wave vector; the wave speed is then c = w/|k|. 

Even simpler to discuss is the first order wave equation 

u t + cu x = 0, (7.1) 

which is trivially solved by characteristics to give 

u = f(x-ct), (7.2) 

representing a wave of speed c. The idea of finding characteristics generalises to 
systems of the form 

Au t + Bu x = 0, (7.3) 

where u € R n and A and B are constant n x n matrices. We can solve this system 
as follows. The eigenvalue problem 

A,4w = Bw (7.4) 

will in general have n solution pairs (w, A), where each value of A is one of the roots 
of the n-th order polynomial 

det(XA -B) = 0. (7.5) 
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Suppose the n values of A, A;, are distinct (which is the general case); then the 
corresponding W; are independent, and the matrix P formed by the eigenvectors as 
columns (i.e., P = (wi, w„)) satisfies BP = APD, where D is the diagonal matrix 
diag(Ai, X n ). P is invertible, and if we write v = P _1 u, then APv t + BPv x = 0, 
whence v 4 + Dv x = 0, and the general solution is 



u 



(7.6) 



representing the superposition of n travelling waves with speeds A^. This procedure 
works providing A is invertible, and also if all the A, are real, in which case we say 
the system is hyperbolic. 

More precisely, we can use the above prescription to solve the nonlinear equation 

Au t + Bu x = c(x,t,u), (7.7) 

where we allow A and B to depend on x and t also. The diagonalisation procedure 
works exactly as before, leading to 

A-(Pv) + B — (Pv) = c[x,t,Pv]; (7.8) 

now, however, A, w and therefore also P will depend on x and t. Thus we find 

v t + Dv x = P~ l A~ l c - [P~ l P t + DP- x P x \v, (7.9) 

and the components of v can be solved as a set of coupled ordinary differential 
equations along the characteristics dx/dt = A,. 

If A and B depend also on u, the procedure is less clear for systems. However, 
the method of characteristics always works in one dimension, so we now return our 
attention to this case. Consider as an example the nonlinear evolution equation 

u t + uu x = 0, (7-10) 

to be solved on the whole real axis. The method of characteristics leads to the 
implicitly defined general solution 

u = f(x-ut), (7.11) 

which is analogous to (7.2), and represents a wave whose speed depends on its am- 
plitude. Thus higher orders of u propagate more rapidly, and this leads to the wave 
steepening depicted in figure 7.1. 

In fact, it can be seen that eventually u becomes multi- valued, and this signifies a 
break down of the solution. The usual way in which this multi- valuedness is avoided 
is to allow the formation of a shock, which is a point of discontinuity of u. the 
characteristic solution applies in front of and behind the shock, and the characteristics 
intersect at the shock, whose propagation forwards is described by an appropriate 
jump condition: see figure 7.2. 



93 



u 




X 

Figure 7.1: Nonlinearity causes wave steepening. 



▲ u 



► 




► 

X 

Figure 7.2: Intersection of characteristics leads to shock formation. 
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This seemingly arbitrary escape route is motivated by the fact that evolution 
equations such as (7.10) are generally derived from a conservation law of the form 



f B pdx = -[q] B A , (7.12) 

J A 



d L r B 
dt 

where the square- bracketed term represents the jump in q between A and B; (7.10) 

2 1 



might be derived from this equation with p = u and q = \u 2 , for example. The 



deduction of the point form 

1+fH ™ 

from (7.12) requires the additional assumption that p and q are continuously difleren- 
tiable; however, it is possible to satisfy (7.12) at a point of discontinuity of p. Suppose 
p is discontinuous at x = xs(t), and denote the jump in any quantity r across the 
shock as [r]t = r(xs+,t) — r(xs—,t). Then let A = xs(t)— and B = x${t + 5t)+. 

rB 

From first principles we have the change of / pdx in the time interval 5t as 

J A 

fB 

5 / pdx pa [p] + 5x s pa [q] + St, (7.14) 

J A 



and from this it follows that 

[p] 

For the particular case of (7.10), we then find 



is = (7.15) 



[^] + 

[«]: 



x s = [ -h^ = \(u + + u_). (7.16) 



An example 

We illustrate how to solve a problem of this type by considering the initial value 

u = u (x) = - 1 , at t = 0. (7.17) 
1 + x 2 

The implicitly defined solution is then 
or, in characteristic form, 

« = «o(0=j^2» x = £ + ut. (7.19) 

This defines a single- valued function so long as u x is finite. Differentiating (7.19) 
leads to 

l + tu' (0' 



Ux = q 7T7TTt\i I 7 - 20 ) 
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X 



Figure 7.3: Characteristic diagram indicating shock formation. 



and this shows that u x — > — oo as t — > t c = min [ r-rrr]- Since — u' Q = 2£/(l + £ 2 ) 2 , 

Z-.u' <o u (£) 

we find the relevant value of £ is l/y/3, and thus t c = 8/ 3 a/3 and the corresponding 
value of x is x c = a/3- Thus (7.18) applies while t < t c = 8/3a/3, and thereafter the 
solution also applies in x < xs(t) and x > xs(t), where 

x s = l[u(xs+) + u(x s -)], (7.21) 

with 

x s = \/3 at t = (7.22) 

As indicated in figure 7.3, the characteristics intersect at the shock, and it is geomet- 
rically clear from figure 7.1, for example, that u + and u_ are the largest and smallest 
roots of the cubic (7.18). An explicit solution for xs is not readily available, but it is 
of interest to establish the long term behaviour, and for this we need approximations 
to the roots of (7.17) when t » 1. 
We write the cubic in the form 

x , 1/1- u\ 1/2 

"-tM— ) • (7 ' 23) 

We know that u < 1, and we expect the largest root, at least, to tend to infinity as 
t — > oo. In that case u ~ x/t if u = 0(1), and the next corrective term gives 

x , i ft - xy/ 2 , mnA . 

• <7 - 24) 
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Figure 7.4: Determination of W(X). 

This evidently gives the upper two roots for x < t (since they coalesce at u = 1 when 
x — t). For large x, the other root must have a<l, and in fact 

(7.25) 

x z 

in order that (7.23) not imply (7.24). Alternatively, (7.25) follows from consideration 
of (7.18) in the form 

tV - 2xtu 2 + {x 2 + l)u - 1 = 0, (7.26) 

providing x ^> t 1 / 3 . 

(7.25) gives the right hand nose of the (rightmost) curve in figure 7.1. The left 
nose can be determined from the observation that the approximation that u ~ x/t 
breaks down (from (7.24)) when x ~ t 1 / 3 , which is also where (7.25) becomes invalid. 
This suggests writing 

u= X t W{X\ X = £ j , (7.27) 
and then W(X) is given approximately, for large t, by 

W{W-l) 2 = ^, (7.28) 

and for X = 0(1) there are three roots providing X > 3/2 2/3 ; at X = 3/2 2/3 , the two 
lower roots coalesce at W — |: this is the left nose of the curve. 

As X becomes large, the upper two roots approach W = 1, x/t, while 

the lower approaches zero, specifically W & 1/X 3 , l/x 2 : see figure 7.4. 

Thus these roots match to the approximations when x ~ t. As X becomes small, the 
remaining root is given by W ~ 1/X, i.e. u ~ l/t 2//3 , and (7.26) shows that this is 
the correct approximation as long as |x| <C t 1 / 3 . The situation is shown in figure 7.5). 
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L u 




x=3t 1/3 /2 2/3 



Figure 7.5: Large time solution of the characteristic solution. 

In order to determine the shock location xs, we make the ansatz that t 1//3 Cxj < 
t, i. e., that the shock is far from both noses. In that case 



u + « u- « — , (7.29) 



and at leading order we have 



whence 



x s « (7.30) 



~ at 1/2 , (7.31) 



confirming our assumption that t 1//3 <C C t 

To determine the coefficient a, we may use the equal area rule, which follows from 
conservation of mass, and implies that the two shaded areas in figure 7.5 are equal. 
We use (7.27) for the left hand area, and (7.24) for the right hand area, then 

/•at 1 / 2 x r t o ft — r\ l l 2 

jUvi^W-"- WV*»Li{— ) *• (7 ' 32) 

where W + and W- are the middle and lowest roots of (7.28), as shown in figure 7.4. 
We write x = t 1//2 £ in the left integral and x = tr) in the right, and hence we deduce 
that 

1/2 

d V = 7T. (7.33) 
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7.1 Traffic flow 



Although the presence of a shock for (7.10) is entirely consistent with the derivation 
of the equation from an integral conservation law, nature appears generally to avoid 
discontinuities and singularities, and it is usually the case that in writing (7.10), we 
have neglected some term which acts to smooth the shock, so that the change of u is 
rapid but not abrupt. 

A simple example of this arises in the theory of traffic flow on a road. On a 
one-lane carriageway, the density of cars is p(x,t) (number per unit length, idealised 
as a continuum, which should be a reasonable assumption for road lengths of many 
cars); x is distance along the carriageway, and t is time. If v(x,t) is the local traffic 
vehicular speed, then a conservation law for cars is 

| + |W = 0, (T.34) 

assuming no cars leave or join the carriageway. The car speed must be prescribed, and 
a simple and sensible prescription is to take v = v(p): a driver drives at a speed which 
depends on density. More specifically, we might suppose that there is a maximum 
speed of vq when there is no traffic (p — 0), but the speed decreases as p increases, 
reaching zero when the distance between cars is zero, corresponding to a maximum 
density p m say. A simple recipe which satisfies these constraints is 

(7.35) 

By scaling the variables suitably, we may take vo = 1, p m = 1, and then we have 

Pt + c{p)p x = 0, (7.36) 

where 

c{p) = l-2p. (7.37) 

This is equivalent to (7.10), with u — 1 — 2p, and therefore we can in general expect 
shocks to form. Since < p < 1, a difference here is that we can have c < for p > |. 

Suppose first that p < \. A local reduction in p is like a bump for u = 1 — 2p, 
and will form a forward propagating shock where p + > p_; cars driving through the 
low density precursor will suddenly hit a traffic jam at high density. Equivalently, a 
local rise in density leads to a forward-propagating shock with p + > p_. 

On the other hand, suppose that p > \. With w = 2p — 1, then w t — ww x = 0, 
and a similar discussion arises, with the direction of x reversed. Thus a local rise in 
p causes a shock to form in which, again, p + > p_, and the shock now propagates 
backwards. In all cases, the individual driver experiences the shock as a sudden 
reduction of speed (the car speed is 1 — p and is larger than the wave speed 1 — 2p: 
thus cars approach and travel through the shock). 

The occurrence of shocks is a common experience to anyone who has driven on a 
motorway, and can be caused by lane closures, motorway junctions, or speed restric- 
tions. Figure 7.6 shows a beautiful series of shock waves propagating backwards in 




99 



morning rush hour traffic on the M25. Beautiful, unless you happen to be in one of 
the vehicles. The diagram actually suggests that the waves arise through an instabil- 
ity, much like roll waves on rivers, but these cannot be described within the confines 
of the simple kinematic wave model (7.36). 

In general (but not quite always), collisions do not occur, and one may ascribe 
this to a smoothing mechanism; the question arises, what could this be? 

Drivers arguably adjust their speed not only according to the inter-car distance 
(the direct density effect), but also depending on how much traffic they see ahead. A 
common exhibition of this is the tailgater: driving at the speed limit on a motorway 
in the outside lane will attract tailgaters if there is little traffic ahead, but not (as 
much) if the road is full. That is, at a constant v, p may increase above its ideal value 
if p in front is less, i.e. if dp/dx < 0. A simple modification for v which describes 
this is (in scaled units) 

v = 1- p - np x , (7.38) 
where k is constant. From this we can deduce the advection-diffusion equation 

£ + <>-*>£- 4 ('£)• ^ 

and the non-local dependence of v on p x is represented as a (nonlinear) diffusion term. 

What is the effect of this on the structure of the solutions? If k is small, we should 
suppose that it is not much, so that shocks would start to form. However, the neglect 
of the diffusion term becomes invalid when the derivatives of p become large. In fact, 
the diffusion term is trying to do the opposite of the advective term. The latter is 
trying to fold the initial profile together like an accordion, while the former is trying 
to spread everything apart. We might guess that a balanced position is possible, in 
which the nonlinear advective term keeps the profile steep, but the diffusion prevents 
it actually folding over (and hence causing a discontinuity), and this will turn out to 
be the case. 

If we suppose instead of (7.38) that 

v = l-p--p x , (7.40) 
P 

which might represent the fact that the strength of tailgating becomes more severe, 
the emptier the road ahead, then the diffusion term becomes linear, and we have 
Burgers' equation for u = 1 — 2p: 

u t + uu x = Ku xx . (7.41) 

7.2 Burgers' equation 

Shock structure 

We suppose k<1, and that 0, and a shock forms at x — xs(t). Our aim 

is to show that (7.41) supports a shock structure, i.e. a region of radical change near 
xs for U- to u+. 



100 




Figure 7.6: Rush hour traffic waves propagating backward on the M25. Horizontal 
axis is minutes from midnight, vertical axis is "outstation number". Separation of 
outstations is 500 m, so the diagram shows 17.5 km of motorway between about 6.40 
a. m. and a little after 11 a. m.: this is a typical weekday picture. The picture is built 
from 1 minute averages of the traffic speed (in lane 2) on the clockwise carriageway 
of the M25 in the neighbourhood of the M3 and M4 interchanges. (The seed point of 
the stop and go waves near the top of the picture is the M4 interchange.) Colour is 
speed in km h _1 , according to colour bar next to picture. Black lines are "simulated 
vehicle trajectories", meaning their gradient at any point is the 1 minute average 
vehicle speed. The picture was made by Eddie Wilson in Matlab from data belonging 
to the Highways Agency. The data is captured by their "MIDAS system" which is 
part of the "Controlled Motorways" project. 
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To focus on the shock, we need to rescale x near xs, and we do this by writing 

x = x s (t) + kX. (7.42) 

We use the chain rule to rewrite the equations in terms of t and X derivatives. To 
be precise, we consider the change of variables t,x — > r,X, where (7.42) defines X 
in terms of t and x, and r = t. Elementary application of the chain rule for partial 
derivatives then shows that 

d^ = d^_x l J)_ d = 1 d 

dt dr k dX' dx kOX' { ' ' 

Burgers' equation becomes 

nu t - x s u x + uu x = u X x, (7.44) 

where we replace r by t since these are in fact the same. We expect the characteristic 
solution (with k = 0) to be approximately valid far from x s , and so appropriate 
conditions (technically, these are matching conditions) are 

u — y u± as X ->• ±oo, (7-45) 

and we take these values as prescribed from the outer solution (i.e., the solution of 

u t + uu x = 0). 

Since k«1, (7.44) suggests that u relaxes rapidly (on a time scale t ~ k « 1) 
to a quasi-steady state (quasi-steady, because u + and m_ will vary with t) in which 

-x s u x + uu x « u X x, (7.46) 

whence 

K - x s u + \u 2 Pa u x , (7.47) 
and prescription of the boundary conditions implies 

K = x s u + - \u\ = x s u- - \u 2 _, (7.48) 

whence 

\-u 2 ] + 

is = (7.49) 

which is precisely the jump condition we obtained in (7.16). The solution for u is 
then 

u = c - (u- - c) tanh \{u- - c)x] , (7.50) 

where c = x$- 



102 



7.3 The Fisher equation 



In Burgers' equation, a wave arises as a balance between nonlinear advection and 
diffusion. In Fisher's equation, 

u t = it(l - it) + u xx , (7-51) 

a wave arises as a mechanism of transferring a variable from an unstable steady state 
(u = 0) to a stable one (u — 1). Whereas Burgers' equation balances two transport 
terms, Fisher's equation balances diffusive transport with an algebraic source term. 
It originally arose as a model for the dispersal of an advantageous gene within a 
population, and has taken a plenary role as a pedagogical example in mathematical 
biology of how reaction (source terms) and diffusion can combine to produce travelling 
waves. 

We pose (7.51) with boundary conditions 

u — y 1, x — > — oo, 

u 0, x ->■ +oo. (7.52) 

It is found (and can be proved) that any initial condition leads to a solution which 
evolves into a travelling wave of the form 

u = f(0, Z = x-ct, (7.53) 

where 

/" + c/' + /(l-/) = 0, (7.54) 

and 

/(oo) = 0, /(-oo) = l. (7.55) 
In the (/, g) phase plane, where g = — /', we have 

f = -9, 

9' = f(l-f)-cg, (7.56) 

and a travelling wave corresponds to a trajectory which moves from (1,0) to (0,0). 
Linearisation of (7.56) near the fixed point (/*, 0) via / = /* + F leads to 




(7.57) 



with solutions e A ^, where A 2 + cA + (1 - 2/*) = 0. We anticipate c > 0; then (1,0) 
is a saddle point, while (0,0) is a stable node if c > 2 (and a spiral if c < 2). For 
c > 2, a connecting trajectory exists as shown in figure 7.7: in practice the minimum 
wave speed c = 2 is selected. (Connecting trajectories also exist if c < 2, but because 
(0,0) is a spiral, these have oscillating tails as u — > 0, which are unstable and also (for 
example, if u is a population) unphysical.) 
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Figure 7.7: Phase portrait of Fisher equation, (7.56), for c = 2. Note how close 
the connecting trajectory (thick line) is to the g nullcline. This is why the large c 
approximation is accurate for this trajectory. 

Explicit solutions for (7.54) are not available, but an excellent approximation is 
easily available. We put 

£ = cE, (7.58) 

so 

*/" + /' + /(I -/) = 0, (7.59) 
with v = 1/c 2 = 1/4 for c = 2. Taking i/Cl and writing / = f + vfi + . . ., we have 

/o + /o(l"/o) = 0, 
fi + (1 - 2/ )/i = -tf, (7.60) 



and thus 



Also, noting that 1 - 2f = -fX/ft (differentiate (7.60)i), 

/i=/o(l-/o)ln[/o(l-/o)], (7.62) 
and so on. Even the first term gives a good approximation, and even for c = 2. 

7.3.1 Stability 
7.4 Solitons 

The Fisher wave is an example of a solitary travelling wave. Another type of solitary 
wave is the soliton, as exemplified by solutions of the Korteweg-de Vries equation 

u t + uu x + u xxx = 0. (7.63) 
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This has travelling wave solutions u = /(£), £ = x — ct, where 

f" + //' - cf = 0, (7.64) 
and solitary waves with / — > at ±oo satisfy the first integral 

f" + lf 2 -cf = 0, (7.65) 

and thus 

\f 2 + If - \cP = 0, (7.66) 

with solution 

/=§csech 2 (^). (7.67) 

Thus there is a one parameter family of these solitary waves, and they are called 
solitons, because they have the remarkable particle-like ability to 'pass through' each 
other without damage, except for a change of relative phase. Despite the nonlinear- 
ity, they obey a kind of superposition principle. Soliton equations (of which there 
are many) have many other remarkable properties, beyond the scope of the present 
discussion. 

Some understanding of the solitary wave arises through an understanding of the 
balance between nonlinearity (uu x ) and dispersion (u xxx ). The dispersive part of 
the equation, u t + u xxx = 0, is so called because waves exp[ik(x — ct)] have wave 
speed c = —k 2 which depends on wavenumber k; waves of different wavelengths 
(2ir/k) move at different speeds and thus disperse. On the other hand, the nonlinear 
advection equation u t +uu x has a focussing effect, which (from a spectral point of view) 
concentrates high wave numbers near shocks (rapid change means large derivatives 
means high wavenumber) . So the nonlinearity tries to move high wavenumber modes 
in from the left, while the dispersion tries to move them to the left: again a balance 
is struck, and a travelling wave is the result. 

7.5 Snow melting 

An example of some of the ideas presented so far occurs in the study of melting 
snow. In particular, it is found that pollutants which may be uniformly distributed 
in snow (e. g. SO2 from sulphur emissions via acid rain) can be concentrated in melt 
water run-off, with a consequent enhanced detrimental effect on stream pollution. 
The question then arises, why this should be so. 

The model we use is based on the principles of groundwater flow. Suppose we have 
a snow pack of depth h in < z < h, where z is a coordinate pointing downwards 
from the snow surface. Snow is a porous aggregate of ice crystals, and meltwater 
formed at the surface can percolate through the snow pack to the base, where run-off 
occurs. (We ignore effects of re- freezing of meltwater) . The flux of water u downwards 
is given by Darcy's law 

(7.68) 



u 



k 



dp 

dz 



+ 99 
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where u is measured as a velocity, and represents volume per unit area per unit time; 
p is the (pore) water pressure, p is density, g is gravity, p is viscosity, k is permeability. 
The permeability k is related to the saturation S, and we will assume 

k = k S 3 , (7.69) 

where the saturation S is the volume fraction of the pore space which is occupied by 
water. For 5=1, the snow is fully saturated (no air), for S = it is fully dry (no 
water). 

Conservation of the liquid water implies 

4MH> «■«> 

where cf) is the porosity (volume of pore space per unit volume of snow). Finally, the 
water pressure is related to the air pressure (p a , taken as constant) by the capillary 
pressure 

Pc=Pa~ P, (7.71) 

and this is a function of S: we take 

Pc(S) = H> Q - S) , ( 7 - 72 ) 

based on typical experimental results. 

Suitable boundary conditions in a melting event might be to prescribe the melt 
flux at the surface: 

u = uq at z = 0. (7.73) 

If the base is impermeable, then 

u = at z = h. (7.74) 

This is certainly not realistic if S reaches 1 at the base, since then ponding must 
occur and presumably melt drainage via a channelised flow, but we examine the 
initial stages of the flow using (7.74). Finally, we suppose S = at t = 0. Again, this 
is not realistic in the model (it implies p c = oo) but it is a feasible approximation to 
make. 

Simplification of this model now leads to the Darcy- Richards equation in the form 



dS 3k pg 2 dS k p d 
0^— H b — = 



dt p dz /j, dz 



(7.75) 



which, as we see, is a convective-diffusion equation of Burgers type. The quantity 
K = kopg/p is known as the saturated hydraulic conductivity; it is a velocity, and 
represents the highest rate at which water can flow through the snow steadily under 
gravity. 

To scale the equation, we note that S is dimensionless and of 0(1) by definition. 
We choose scales for z and t: 

z~h, t~^; (7.76) 
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the model then becomes 

OS 



55 352 as _ d_ 

dt dz dz 



(7.77) 



where 

K = (7.78) 
pgh 

if we choose po — 1 kPa, p = 10 3 kg m -3 , g = 10 m s -2 , /i = 1 m, then /c = 0.1. 
it follows that (7.77) has a propensity to form shocks, these being diffused by the 
term in k over a distance O(k) (by analogy with the shock structure for the Burgers 
equation). 

We want to solve (7.77) with the initial condition 

S = at t = 0, (7.79) 

and the boundary conditions 

S 3 - kS(1 + S 2 )^ = ^ on z = 0, (7.80) 
oz K 

and 

5 3 -/c5(l + 5 2 )— = at z=l. (7.81) 

(JZ 



Roughly, for « <C 1, these are 



S = So at z = 0, 
5 = at z = l, (7.82) 



where = (uo/K) 1 ^ 3 , which we initially take to be 0(1) (and < 1, so that surface 
ponding does not occur). 

Neglecting k, the solution is the step function 

S — Sq, z < Zf, 

S = 0, z > z f , (7.83) 
and the shock front at Zf advances at a rate Zf, given from the jump condition 

*/ = fj£ = S ^ < 7 - 84 ) 

In dimensional terms, the shock front moves at speed uo/4>So, which is in fact obvious 
(given that it has constant S behind it). 

The shock structure is similar to that of Burgers' equation. We put 

z = z f + kZ, (7.85) 

and S rapidly approaches the quasi-steady solution S(Z) of 

-VS' + 3S 2 S' = [5(1 + S 2 )S']', (7.86) 



107 




Figure 7.8: S(Z) given by (7.90); the shock front terminates at the origin, 
where V — Zf\ hence 

S{1 + S 2 )S' = -S(S 2 -S 2 ), 
in order that S — > So as Z — > — oo, and where we have chosen 

V = sl 

(as S + — 0), thus reproducing (7.84). The solution is a quadrature, 



/ 



s (1 + S 2 )dS 
(S 2 - 52) 



with an arbitrary added constant (amounting to an origin shift for Z). Hence 



(1 + So 2 ) 
2S n 



In 



Sn 



S -S 



= Z. 



(7.87) 
(7.88) 

(7.89) 
(7.90) 



The shock structure is shown in figure 7.8, and one particular feature is notewor- 
thy; the profile terminates where S = at Z = 0. In fact, (7.87) implies that 5 = 
or (7.90) applies. Thus when S given by (7.90) reaches zero, the solution switches to 
5 = 0. The fact that dS/dZ is discontinuous is not a problem because the diffusivity 
5(1 + S 2 ) goes to zero when 5 = 0. In fact, this degeneracy of the equation is a 
signpost for fronts with discontinuous derivatives, and we shall encounter this situa- 
tion again when we study non-linear diffusion. Essentially, the profile can maintain 
discontinuous gradients at 5 = because the diffusivity is zero there, and there is no 
mechanism to smooth the jump away. 

Suppose now that k = 10~ 9 m 2 (a plausible value) and \ij p = 10~ 6 m 2 s _1 , then 
the saturated hydraulic conductivity K = k pg/fi = 10~ 2 m s _1 . On the other hand, 
if a metre thick snow pack melts in ten days, this implies uq ~ 10~ 6 m s _1 . Thus 
Sq = uq/K ~ 10~ 4 , and the approximation 5 w So looks less realistic. With 



S 3 - kS(1 + S 2 ) 



as 

~d~z 



c3 



(7.91) 
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and So ~ 10 2 and k ~ 10 1 , it seems that one should assume S <C 1. We define 



= (?) 



note that (S$/k) 1/2 ~ 0.03, so that (7.91) becomes 

ds 



l + ^-s 2 

K 



dz 



= 1 on z — 0, 



and we have - 10" 3 , /3 = (S /k) 3 / 2 ~ 0.3. 
We neglect the term in Sq/k, so that 



ps — s— « 1 on 2; 
oz 



and substituting (7.92) into (7.77) leads to 

ds „ n 9 ds d 
1 3/3s 2 



dr 



dz dz 



ds^ 
dz 



(7.92) 



(7.93) 



(7.94) 



(7.95) 



A simpler analytic solution is no longer possible, but the development of the solution 
will be similar. The flux condition (7.94) at z = allows the surface saturation to 
build up gradually, and a shock will only form if /3 ^> 1 (when the preceding solution 
becomes valid). 



7.6 Notes and references 

A fundamental reference on the subject of waves in all kinds of media is the book by 
Whitham (1974). An older book which details waves in fluid is the book by Stoker 
(1957). 

Exercises 

7.1 A simple model for the flow of a mixture of two fluids along a tube (e. g., air 
and water) is 

a t + (av) z = 

-a t + [(1 - a)u] z = 

p g [{av)t + [av 2 ) z ] = -ap z , 

p t [{(l ~ a)u} t + {A(l - a)u 2 } z ] = -(1 - a)p z , 

where p is pressure, u and v are the two fluid velocities, a is the volume fraction 
of the fluid with speed v, p g is its density, and p\ is the density of the other 
fluid. Show that there are two characteristic speeds dz/dt = A, satisfying 

(A - uf = (A - l)[u 2 + 2u(X - u)] - s 2 (X - v) 2 , 
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where 



p g (l - a) 



1/2 



PlOL 



Deduce that the characteristic speeds are real if, when Di — 1<1, sCl, 

2 

I .Sl 7/. — 11 ) I 

A > 1 + 



s{u — v) 



u 



In particular, show that the roots are complex if Di — 1 and u ^ v. What does 
this suggest concerning the well-posedness of the model? 

7.2 The function u(x,t) satisfies 

u t + uu x = a(l — u 2 ) 

for — oo < x < oo, where alpha > 0, and with m = uo(x) at t = 0, and 
< uq < 1 everywhere. Show that the characteristic solution can be written 
parametrically in the form 

u (s) + tanhat . . . . 

M = , exp \a(x — s) \ — cosh at + unis) smhat. 

1 + u (s)tanhat' FL v n w 

Sketch the form of the characteristics for an initial function such as uo(s) = 
a/(l + s 2 ). Show that, in terms of s and t, u x is given by 



it. 



[asech 2 at]u' (s) 



[1 + wq(s) tanh at] [a + {«o(s) + auo(s)} tanhat] 



and deduce that a shock will form if u' + a(l + uo) becomes negative for some 
s. Show that if uq = a/(l + s 2 ) and a is small, this occurs if 



a < 



3a\/3 
8 



7.3 Discuss the formation of shocks and the resulting shock structure for the equa- 
tion 

U t + U a U x = £[M^M X ] X , 

where a,/3 > 0, and £<1. (Assume u > 0, and u — > at ±oo.) 
Show that the equation 

7J 4 + ttM x = £Mtt xx 

admits a shock structure joining «_ to a lower value u + , but not one in which 

i + 



the wave speed c 



/[u]±. Why should this be so? 



7.4 Use phase plane methods to study the existence of travelling wave solutions to 
the equation 

u t = u p (l - u q ) + [u r u x ] x , 
when (i) p — 1, q — 2, r — 0; (ii) p — 1, q — 1, r — 1. 
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7.5 Two examples of integrable partial differential equations which admit soliton 
solutions are the nonlinear Schrodinger (NLS) equation 

iu t = \u\ 2 u + u xx , 

and the sine-Gordon equation 

u t t - u xx = sinw. 

Show that these equations admit solitary wave solutions (which are in fact 
solitons). 

7.6 In a model of snow melting, it is assumed that the permeability is k = koS a , 
and the capillary suction is p c (S) = po(S~ 13 — S), where a, j3 > 0, and S is 
the saturation. How does the choice of different values of a and j3 affect the 
formation and propagation of shock waves, and their internal structure? 
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Chapter 8 
Similarity solutions 



Another class of particular solutions of partial differential equations is that of sim- 
ilarity solutions. Like travelling wave solutions, similarity solutions are important 
indicators of solution behaviour. They are appropriate when the statement of the 
differential equation and its associated boundary and initial conditions make no ref- 
erence to any intrinsic scale. One way to recognise when this is the case is to see 
whether the equation and its associated initial/boundary conditions are invariant un- 
der a rescaling of the variables. If this is the similarity solution is appropriate, 
and its form can be inferred from the combinations of variables which are invariant 
under the transformation. 



8.1 The heat equation 

The classical situation where a similarity solution is appropriate is the heat equation 
in one space dimension, 

u t = u xx , (8.1) 
which we choose to solve on (—00, 00) subject to the boundary conditions 

u -)• as x -> ±00, (8.2) 

and an initial condition to be described below. 

The heat or diffusion equation represents the diffusion of a quantity whose density 
is u. For example, it may represent heat transfer, where u is temperature, or salt 
diffusion, where u is salt concentration. A standard kind of problem to consider is 
then the release of a concentrated amount at x = at t = 0. We can idealise this by 
supposing that at t — (with u suitably normalised), 

/oo 
u(x)dx = l. (8.3) 
-00 

This apparently contradictory prescription idealises the concept of a very concentrated 
local injection of u. For example, (8.1) with (8.3) could represent the diffusion of sugar 
in hot (one-dimensional) tea from an initially emplaced sugar grain. (8.3) defines the 
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delta function 5(x), an example of a generalised function. One can think of generalised 
functions as being (denned by) the equivalence classes of well-behaved functions u n 
with appropriate limiting behaviour. For example, the delta function is defined by 
the class of well-behaved functions u n for which 

/oo 
u n (x)f(x)dx^f(0) (8.4) 
-oo 

as n — > oo for all well-behaved f(x). As a shorthand, then, 

/oo 
6(x)f(x)dx = f(0) (8.5) 
-oo 

for any /, but the ulterior definition is really in (8.4). In practice, however, we think 
of a delta function as a 'function' of x, zero everywhere except for a (very) sharp spike 
at x = 0. 

A similarity solution is appropriate for this problem because it has no intrinsic 
space or time scales. The time scale is semi-infinite and the space scale is infinite. 
The equation fixes a diffusive time scale if there is a length scale, but there is none. It 
is in this context that one can expect the solution to look the same at different times 
on different scales, and it is this self-similarity that lends the solution type its name. 
In general, as t varies, then the length scale might vary as £(t) and the amplitude of 
the solution u might vary as U (t). That is, if we look at u/U as a function of it 
will look the same for all t. This in turn suggests that the solution takes the form 



u = U(t)f 



x 

W) 



(8.6) 



and this is one of the forms of a similarity solution. 

It is most often the case that U and £ are powers of t, and the exponents are then 
to be chosen so that the problem has such a solution. This is best seen by example. 
If we denote rj = x/£(t), and substitute the form (8.6) into (8.1), (8.2) and (8.3), then 
we find 



/(±oo) = 0, 

m = o, 

fdrj = 1. (8.7) 



Since / is supposed to be a function of r\ only, all t dependence must vanish from 
this set of equations. The integral constraint requires 

U£ = 1 (8.8) 

(the constant can be arbitrarily chosen), whence the differential equation for / be- 
comes 

/" + ^W)' = 0, (8.9) 
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and we choose ££' = 2; together with £(0) = (which follows from the requirement 
that for x ^ 0, u = at t = 0), this implies 

£ = 2Vi, U = ^=, (8.10) 

and the solution for / satisfying the integral constraint is 

/=^e"" 2 . (8.11) 

Finally, we have the fundamental solution of the heat equation, 

«=— ^=exp(-— V (8.12) 

A feature of similarity solutions is the reduction in order, and an associated re- 
duction in the number of boundary conditions and in the number of independent 
variables. Thus, apart from the requirement that the equation admit a similarity 
solution, we need at least two of the initial/boundary conditions to collapse into one. 
In diffusional problems, this is most commonly the initial condition and the condition 
at infinity. This is the case above, where we required u — > both as x — > ±oo and as 
t — > 0. The partial differential equation is third order in the sense that there is one 
time derivative and two space derivatives (thus requiring three conditions), but the 
similarity equation for / is a second order ordinary differential equation. 

There are two other comments to make about the example above. One is that 
the first integral of the similarity equation is a consequence of the conserved quan- 

tity / udx, whose constancy follows from the equation and boundary conditions. 

Joo 

When there is no such conservation law, then generally the equation cannot be solved 
analytically. 

The second comment concerns the usefulness of similarity solutions. They require 
particular combinations of initial and boundary conditions which, one might suppose, 
are not generally applicable. An important observation about more general problems 
is that they will often develop self-similar characteristics in localised regions of space 
or at large times. 

We can illustrate this simply by solving (8.1) and (8.2) with the general initial 
condition 

u = uo(x) at t = 0. (8.13) 

Because the problem is linear, we can use the fundamental solution (8.12) to generate 
a solution. 1 The initial condition can equally be written as 

/oo 
u (05(x-0d£, (8.14) 
-oo 

and this tells us that the solution satisfying (8.13) is 

1 We are simply using the Green's function for the heat equation here. 



114 



Now the idea is this. Suppose m is localised over some region near the origin; 

for example we might have uo = ^; then at large time, we expect the initial 

1 + x 1 

concentration to have diffused a long way, and to have 'forgotten' details of the initial 
condition. It is because of this that a similarity solution becomes relevant. This can 
be seen by approximating the exact solution (8.15) for t ^> 1. We rewrite the solution 

as 



(e - 2fr) 

At 



(8.16) 



For large t and localised uo, the exponential in the integrand becomes approximately 
one for bounded £ (and for large |£| we suppose uo is negligible), and therefore we 
regain (8.12) in the form 



u 



* (£ U . m )exp(4), «>1. (8.17) 



8.1.1 Error function profile 

Let us turn to a problem consisting of the same heat equation, but now to be solved 
on < x < oo, with 

u = at t = 0, 
u — 1 at x = 0, 

u -> as x -)• oo. (8.18) 

There is again a similarity solution, but the x = boundary gives a scale for u, 
and its size is thus fixed. Invariance of the equation again implies that the similarity 
variable can be taken to be x/2y/t, and therefore we seek a similarity solution of the 
form 

« = /fo), ^=^f t - ( 8 - 19 ) 

Substituting this into the heat equation (8.1), together with the initial/boundary 
conditions (8.18), we find that / satisfies 

f" + 2nf = 0, /(0) = 1, /(oo) = 0. (8.20) 

The solution of this is the complementary error function profile 

/foHerfcq, (8.21) 

where the complementary error function is defined by 

2 roc 



erfc 



v = — e~ s ds; (8.22) 

V 7T Jr) 



it is a monotonically decreasing function as shown in figure 8.1, and is a commonly 
occurring function in diffusion problems (the error function is defined by erf 77 = 
1 — erfc 77) . 
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Figure 8.1: The complementary error function. 



8.2 The porous medium equation 

Gas in a porous medium satisfies the mass conservation law 



p t + V. q = 0, 



(8.23) 



where p is density and q is the mass flux. In a porous medium, Darcy's law (which 
we have already met in section 7.5) relates q to pressure gradient, specifically 



pk 
H 



(8.24) 



where k is the permeability, p is the gas viscosity, and p is pressure. The density p 
of a gas is related to its pressure by a constitutive law, which is usually taken to be 
the perfect gas law 

where M is the molecular weight, T is temperature, and R is the universal gas con- 
stant. In isothermal conditions where the temperature is constant, we thus have 
p oc p, and in suitable dimensionless units, we have the porous medium equation 



Pt = V. [pVp]. 



(8.26) 



This is appropriate for slow gas seepage. 

If the gas flow is rapid, then the temperature alters as the pressure changes, and 
in adiabatic conditions where there is no heat loss, p oc p 7 , where 7 > 1 is the ratio 
of the specific heats; for air, 7 = 1.4. In this case, the appropriate version of (8.26) is 



(H = V. [p m Vpl 



(8.27) 
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where m — 7 > 1. Thus both cases are described by (8.27), differing only in the 
choice of exponent m. In both cases, m > 1. 

Let us consider the seepage of gas in a one-dimensional column, and let us write 
u for p, so that 

u t = {u m u x ) x , (8.28) 
and let us consider the spread of an initially concentrated dose at x = 0, so that 

u — > as x — > ±00, 

u = at t = 0, x ^ 0, (8.29) 

and 



/oo 
m dx = 1 (8.30) 
-00 



for all t. 

We seek a similarity solution of the form 

u = U(t)f(r 1 ), V=^- (8-31) 

The integral constraint requires U = l/£, and then we find 

-r +1 £V)' = [r/T, (8.32) 

where £' = but /' = df jdr\. The initial/boundary conditions become 

/(±oo) = 0, (8.33) 
and the normalisation condition (8.30) is 

/oo 

fd V = l. (8.34) 
-00 

The solution requires to be constant, and it is algebraically convenient to 

choose £ m+1 £' = 2/m, thus 



77 = x 



m 



2(m + 2)t 



ro+2 

(8.35) 



Because u is conserved with these boundary conditions, a first integral of (8.32) exists, 
and is 

rr + -vf = 0, (8.36) 

m 

with the constant of integration being zero (because / — > as 77 — > ±00). Thus either 
/ = 0, or 

/ = H ~ V 2 ] 1/m , (8-37) 
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so that the solution has the form of a cap of finite extent, given by (8.37) (for \rj\ < rj , 
and / = for \n\ > rj . The value of t]q is determined from /f^ / dn = 1, and is 



Vo 



if 

Jo 



*/ 2 rn+2 

cos "> 6 d6 



ro+2 



(8.38) 



The finite extent of the profile is due to the degeneracy of the equation when m > 0. 
(The limit m — > regains the Gaussian solution of the heat equation by first putting 
j] = y^m^C, f = F/y/m, and noting that r] & (irm)~ m ^ 2 as m — > (this last 
following by application of Laplace's method 2 to (8.38)).) The graph of f(rj) is shown 
in figure 8.2. 




Figure 8.2: f{rj) given by (8.37). 



8.3 Snow melting revisited 



Let us return to the snow melting model of chapter 7. Suppose that in the model 
described by (7.94) and (7.95), the parameter (3 <C 1; then the saturation profile 
approximately satisfies 



ds 
dz 



d_ 
dz 



dz 

1 on z 
on z 



0, 
1. 



At least for small times, the model admits a similarity solution of the form 

S = T a f( V ), r, = z/rP, 



(8.39) 



(8.40) 



where satisfaction of the equations and boundary conditions requires 2a = (3 and 
2/5 = 1 = a, whence a = 1/3, /3 — 2/3, and / satisfies 



(//')' -|(/ -2r?/') = 0, 



(8.41) 



2 Laplace's method for the asymptotic evaluation of integrals is described in Carrier, Krook and 
Pearson (1966), for example. 
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and the condition at z = becomes 

-ff = l at r/ = 0. (8.42) 

The condition at z — 1 can be satisfied for small enough r, as we shall see, because 
the equation (8.41) is degenerate, and / reaches zero in a finite distance, 770, say, and 
/ = for 77 > 770 . As 77 = 1/r 2 / 3 at z — 1, then this solution will satisfy the no flux 
condition at z = 1 as long as r < Vq 3 ^ 2 , when the advancing front will reach z — 1. 
To see why / behaves in this way, integrate once to find 

f(f + lv) = -l + £fdv. (8.43) 

For small 77, the right hand side is negative, and / is positive (to make physical sense), 
so / decreases (and in fact /' < —\t])- For sufficiently small /(0) = fo, f will reach 
zero at a finite distance 77 = 770, and the solution must terminate. On the other 
hand, for sufficiently large fo, Jq fdrj reaches 1 at 77 = 771 while / is still positive (and 
/' = —§771 there). For 77 > 771, then / remains positive and /' > — 1?7 (/ cannot reach 
zero for 77 > 771 since Jq fdrj > 1 for 77 > 771). Eventually / must have a minimum and 
thereafter increase with 77. This is also unphysical, so we require / to reach zero at 
77 = 770. This will occur for a range of fo, and we have to select fo in order that 

/ fdri = l, (8.44) 
Jo 

which in fact represents global conservation of mass. Figure 8.3 shows the schematic 
form of solution both for j3 ^> 1 and j3 <C 1. Evidently j3 ~ 1 will have a travelling 
front solution between these two end cases. 

8.4 Notes and references 

An excellent introduction to the subject of partial differential equations from the 
viewpoint of the applied mathematician is the book by Carrier and Pearson (1976). 

Exercises 

8.1 The function u(x,t) satisfies the heat equation 

— ^xx 

on < x < 00, together with the initial condition 

u = at t = 0, 

and the boundary condition 

u — > as x — > 00. 
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Figure 8.3: Schematic representation of the evolution of s in (7.95) for both large and 
small /3. 



Find similarity solutions for u when the boundary condition at x = 0, t > has 
each of the following forms: 

(i) u = 1; 

(ii) u x = -1; 
(hi) uu x = -1; 
(iiii) u x = — u' 3 . 

Show in case (iiii) that no similarity solution exists if /3 — 1. What is it about 
this case that prevents a similarity solution existing? 

8.2 Which of the following equations have similarity solutions? In all cases suppose 
that the domain of solution is < x < oo, and that u = at t = and u = 1 at 
x = 0. If there is a similarity solution, solve the similarity equation if you can. 

(i) u t = xu xx ; 

(ii) u t = (xu x ) x ; 
(hi) u t = u\. 

8.3 Write down the equation satisfied by a similarity solution of the form u = t^f(rj), 
7] = x/t a , for the equation 

u t = (u m u x ) x in < x < oo, 

where m > 0, with u m u x = — 1 at x = 0, u — > as x — > oo, « = at t = 0. 
Show that J °° / CZ77 = 1, and show that in fact / reaches zero at a finite value 
770 • Is the requirement that m > necessary? 
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8.4 u satisfies the equation 



u t = [D(u)u x ] x in < x < oo, 

with u = at x — > oo and t = 0. For a general function D (not a power of it), 
for what kind of boundary condition at x = does a similarity solution exist? 
What if, instead, D = D(u x )l Write down suitable equations and boundary 
conditions for the similarity function in each case. 
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Chapter 9 
Nonlinear diffusion 



9.1 The viscous droplet 

An example of where the nonlinear diffusion equation can arise is in the dynamics of 
a drop of viscous fluid on a level surface. If the fluid occupies < z < h(x,y,t) and 
is shallow, then lubrication theory gives the approximation 

„ d 2 u 

Pz = -pg, (9.1) 

in which u = (u, v, 0) is the horizontal component of velocity, and V is the horizontal 
gradient (d/dx,d/dy,0). With p = at z = h, we have the hydrostatic pressure 
p = pg(h — z), so that = pgV/i, and three vertical integrations of (9. l)i (with 
zero shear stress du/dz = at z = h and no slip u = at z = 0) yield the horizontal 
fluid flux 

q = ^ udz = -^h 2 Vh. (9.2) 
Conservation of fluid volume for an incompressible fluid is h t + V.q = 0, and thus 

/», = ^V.[/» 3 V/»], (9.3) 

corresponding to (8.27) (in two space dimensions) with m = 3. A drop of fluid placed 
on a table will spread out at a finite rate. 

That this does not continue indefinitely is due to surface tension. Rather than 
having p = at z = h (where the atmospheric pressure above is taken as zero), the 
effect of surface tension is to prescribe 

p = 27/c, (9.4) 

where 7 is the surface tension, and k is the mean curvature relative to the fluid droplet 
(i.e. ac > if the interface is concave, as illustrated in figure 9.1). The curvature is 
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Figure 9.1: The surface shown has positive curvature when the radius of curvature is 
measured from below the surface; in this case equilibrium requires p > p a . 



defined as 2k = V.n, where n is the unit normal pointing away from the fluid (i.e., 
upwards) . At least this shorthand definition works if we define 



n 



(-h x ,-h y ,l) 
[l + lV/il 2 ] 1 ^' 



thus 



2k 



(9.5) 



(9.6) 



_{i + \vh\ 2 yi\ 

It is less obvious that it will work more generally, since there are many ways of defining 
the interface as cf)(x, y, z) = and thus n = V0/|V</>| (that is (9.5) uses <fi = z — h); 
but in fact it does not matter, since we may generally take cf) = (z — h)P, so that 
V</> = (— h x , —h y , 1) on z = h, and Vc/>/|V</>| is the same expression as in (9.6). 



For shallow flows, we replace p = on z = h by p = — 7V h there, and thus 

p ps pg(h — z) — ^fV 2 h, 



and (via (9.2)), (9.3) is modified to 



h t = V. 



3/i 



V{pgh - 7 V 2 M 



(9.7) 



(9.8) 



The fourth order term is also 'diffusive', insofar as it is a smoothing term: this 
is most easily seen by considering the fate of modes e at+lkx for the linear equation 
h t = —h xxxx : a = —k 4 : high wave number (high gradient) modes are rapidly damped. 
Surface tension can thus also act to smooth out shocks. The effect of surface tension 
relative to the diffusional gravity term is given by the Bond number 



Bo = 



P9t 

7 



(9.9) 
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where I is the lateral length scale of the drop. This is the (only) dimensionless 
parameter which occurs when (9.8) is written dimensionlessly. 

For a drop released on a table, (9.8) still predicts unending dispersal, but if the 
full nonlinear curvature term (9.6) is kept, then a steady state will exist, and surface 
tension keeps the drop of finite extent. 

9.2 Advance and retreat: waiting times 

The similarity solution (8.37) predicts an infinite slope at the margin (where / = 0) 
if m > 1 (and a zero slope if m < 1). If one releases a finite quantity at t — 0, then 
one expects the long time solution to be this similarity solution. The question then 
arises as to how this similarity solution is approached, in particular if the initial drop 
has finite slope at the margin. 

This question can be addressed in a more general way by studying the behaviour 
near the margin x = xs{t) of a solution h(x,t) of the same equation as in (8.28), 

h t = (h m h x ) x . (9.10) 

Suppose that h ~ c(xs — x) v for x near xs- Then satisfaction of (9.10) requires 

x s ~ c m \v{m + 1) - l](x s - xY m -\ (9.11) 

Note that the similarity solution (8.37) has xs finite when v = 1/m, consistent with 
(9.11), and more generally we see that the margin will advance at a rate xs ~ c m /m 
if h ~ c{x s - x) l l m . 

Suppose now that m > 1, and we emplace a drop with finite slope, v — 1. Then 
the right hand side of (9.11) is zero at x = xs, and thus is = 0: the front does not 
move. What happens in this case is that the drop flattens out: there is transport 
of h towards the margin, which steepens the slope at xs until it becomes infinite, at 
which point it will move. This pause while the solution fattens itself prior to margin 
movement is called a waiting time. 

Conversely, if m < 1, then the front moves (forward) if the slope is zero, v = 1/m. 
If the slope is finite, v — 1, then (9.10) would imply infinite speed. An initial drop of 
finite margin slope will instantly develop zero front slope as the margin advances. 

(9.11) does not allow the possibility of retreat, because it describes a purely dif- 
fusive process. The possibility of both advance and retreat is afforded by a model of 
a viscous drop with accretion, one example of which is the mathematical model of an 
ice sheet. Essentially, an ice sheet, such as that covering Antarctica or Greenland, 
can be thought of as a (large) viscous drop which is nourished by an accumulation 
rate (of ice formed from snow) . A general model for such a nourished drop is 

ht = (h m h x ) x + a, (9.12) 

where a represents the accumulation rate. Unlike the pure diffusion process, (9.12) 
has a steady state 

(x 2 -x 2 )^, (9.13) 



h 



(m + l)a 
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where x must be prescribed. (In the case of an ice sheet, we might take x to be at 
the continental margin.) (9.13) is slightly artificial, as it requires a = for x > x , 
and allows a finite flux —h m h x = axo where h = 0. More generally, we might allow 
accumulation and ablation (snowfall and melting), and thus a = a(x), with a < for 
large In that case the steady state is 



h 

where the balance is 



(m + 1) / sdx 

J X 



m+1 

(9.14) 



/ at 

Jo 



,dx, (9.15) 
and Xq is defined to be where accumulation balances ablation, 

adx = 0. (9.16) 



/* 

Jo 



This steady state is actually stable, and both advance and retreat can occur. 
Suppose the margin is at xs, where a = a$ = —\a>s\ ( a s < 0, representing ablation). 
If we put h « c(xs — x) u , then (9.12) implies 

vcx s (x s - xy- 1 ~ vc m+1 [v{m + 1) - l](x s - x )M m +V-% - \a s \, (9.17) 

and there are three possible balances of leading order terms. 
The first is as before, 

x s ~ c m [v{m + 1) - l](x s - xY m -\ (9.18) 

and applies generally if v < 1. Supposing m > 1, then we have advance, xs ~ c m /m 
if v = 1/m, but if v > 1/m, this cannot occur, and the margin is stationary if 
1/m < v < 1. If v — 1, then v{m + 1) — 2 = m — 1 > 0, so that 

x S ™-\as\/c, (9.19) 

and the margin retreats; if v > 1, then instantaneous adjustment to finite slope and 
retreat occurs. 

The ice sheet exhibits the same sort of waiting time behaviour as the viscous drop 
without accretion. For 1/m < v < 1, the margin is stationary, and if x$ < xq then 
the margin slope will steepen until v = 1/m, and advance occurs. On the other hand, 
if xs > xo, then the slope will decrease until v — 1, and retreat occurs. In the steady 
state, a balance is achieved (from (9.14)) when v = 2/(m + 1). 



9.3 Blow-up 

Further intriguing possibilities arise when the source term is nonlinear. An example 
is afforded by the nonlinear (diffusion) equation 

u t = u xx + Xe u , (9.20) 
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Figure 9.2: Maximum values of u, u(0), as a function of the parameter A. Blow-up 
occurs if A > 0.878. 



which arises in the theory of combustion. Indeed, as we saw earlier, combustion occurs 
through the fact that multiple steady states can exist for a model such as (7.74), and 
the same is true for (9.20), which can have two steady solutions. In fact, if we solve 



u 



Xe u = with boundary conditions a = on i = ±1, then the solutions are 



u = 2 In 



ylsech 




(9.21) 



where A = exp[u(0)/2], and A satisfies 

A — cosh 




(9.22) 



which has two solutions if A < 0.878, and none if A > 0.878: the situation is depicted 
in figure 9.2. If we replace e u by exp[«/(l + eu)], e > 0, we regain the top (hot) 
branch also. 

One wonders what the absence of a steady state for (9.20) if A > A c implies. 
The time-dependent problem certainly has a solution, and an idea of its behaviour 
can be deduced from the spatially independent problem, u t = Xe u , with solution 
u = ln[l/{A(to _ t)}]\ u reaches infinity in a finite time. This phenomenon is known 
as thermal runaway, and more generally the creation of a singularity of the solution 
in finite time is called blow-up. Numerical solutions of the equation (9.20) including 
the diffusion term show that blow-up still occurs, but at an isolated point; figure 9.3 
shows the approach to blow-up as t approaches a critical blow-up time t c . 

In fact, one can prove generally that no steady solutions exist for A greater than 
some critical value, and also that in that case, blow-up will occur in finite time. To 
do this, we use some pretty ideas of higher grade mathematics. 
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3.56350943288956 




3.56383223992489 


- 


3.56383806154693 



















-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 
X 

Figure 9.3: Solution of u t = u xx + e u on [—1, 1], with u = at x = — 1, 1 and t = 0. 
The solution is shown for four times close to the blow-up time t c = 3.56384027594971. 
The many decimal places should not be treated too seriously, but they do indicate 
the logarithmic suddenness of the runaway. 

Suppose we want to solve the more general problem 

ut = V 2 U + Xe u in ft, (9.23) 

with u = in the boundary dQ, and u = at t = (these conditions are for 
convenience rather than necessity). We will be able to prove results for (9.23) which 
are comparable to those for the ordinary differential equation version (cf. (7.81)) 

w = -frw + Xe w , (9.24) 

because, in some loose sense, the Laplacian operator V 2 resembles a loss term. 

More specifically, we recall some pertinent facts about the (Helmholtz) eigenvalue 
problem 

V 2 (p + (Mp = in ft, (9.25) 

with (j) = on dfl. There exists a denumerable sequence of real eigenvalues < \i\ < 
fx 2 ■ ■ with /j, n — > oo as n — > oo, and corresponding eigenfunctions </>i,</>2, ■■■ which 
form an orthonormal set (using the I? norm), thus 

</>j) = f (PrfjdV = 5ij, (9.26) 

where Sij is the Kronecker delta (= 1 if i — j, if i ^ j). These eigenvalues satisfy a 
variational principle of the form 

m = min I IV0IW, (9.27) 
Jn 
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where (j) ranges over functions of unit norm, ||</>|| 2 = {/ ^dV} 1 ^ 2 = 1, which are 
orthogonal to <f>j for j < i; (more generally fa = min {/ |V</>| 2 / / </> 2 } if </> is not 
normalised on to the unit sphere ||</>|| 2 = 1). In particular 



lin f |V</>|W, (9.28) 
h=iJn 



Ui = mm 

and the corresponding </>i is of one sign, let us say positive. 

We take the inner product of the equation (9.23) with cf>i and divide by / 0i dV; 
defining 

v{t) = j^dV=L ud ^ (9 - 29) 

where dec = <pi dV j Jq4>i dV is a measure on fl (with J n dco = 1), and using Green's 
theorem, we find 

v = A f e u du- fav, (9.30) 
Jn 

and the equation for v is close to the ordinary differential equation (9.24). 

Now we use Jensen's inequality. This says that if we have an integrable function 
g on Q and a convex function f(x) (i.e. one that bends upwards, /" > 0), then 



/ 



/ gdJ< [ f[g]du (9.31) 
Jn J Jn 



for any measure u on Q such that f n clou = 1. We have chosen u to be so normalised, 
and e u is convex: thus 

J exp(u) doj > exp udoj — e v , (9.32) 

so that 

v > Xe v - p x v. (9.33) 

It is now easy to prove non-existence and blow-up for A greater than some critical 
value A c . Firstly, u must be positive, and hence also v. (For suppose u < 0: since 
u = at t = and on dQ, then u attains its minimum in Q at some t > 0, at which 
point u t < 0, u xx > 0, which is impossible, since then u t — u xx = \e w < 0.) For any 
v, e v > ev, thus v > (Ae — n\)v. In a steady state, v = 0, and v > (clearly u = is 
not a solution), so this is impossible if 

A > fa/e. (9.34) 

This implies non-existence of a solution for A > A c , where A c < \x\je. 
In a similar vein, if A > y^i/e, then 

v > iii[e v - 1 - v], (9.35) 

and v > w, where 

w = iiiie™' 1 - w), w(0) = 0. (9.36) 
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(This is a standard comparison argument: v = w at t = 0, and v > w there, so v — w 
is initially positive. It remains so unless at some future time v — w reaches zero again, 
when necessarily v — w < — which is impossible, since v > w whenever v = w.) 
But w — > oo in finite time (w > so that w — > oo as t increases, and as w — > oo, 
e _ ™i£; « yue -1 , so e~ w reaches zero in finite time); therefore also v reaches infinity in 
finite time. Finally 

v— wdw < supw, (9.37) 
Jn n 

since f n du = 1: hence u — >• oo in finite time. 

In fact m — > oo at isolated points, and usually at one isolated point. As blow-up 
is approached, one might suppose that the nature of the solution in the vicinity of 
the blow-up point would become independent of the initial (or boundary) conditions, 
and thus that some form of similarity solution might be appropriate. 

This is indeed the case, although the precise structure is rather complicated. 
We examine blow-up in one spatial dimension, x. As a first guess, the logarithmic 
nature of blow-up in the spatially independent case, together with the usual square- 
root behaviour of the space variable in similarity solutions for the diffusion equation, 
suggests that we define 

T = -ln(t -t), 77= ^1^1/2 , u = -ln[X(t -t)] + g(r 1 ,r), (9.38) 

where blow-up occurs at x = xq at t = to] hence g satisfies 

9t = g m ~ lv9v + e 9 - 1. (9.39) 

The natural candidate for a similarity solution is then a steady solution g(rj) of (9.39), 
satisfying 

(/" - | W ' + (e» - 1) = 0, (9.40) 
and matching to a far field solution u(x,t ) would suggest 

<7~— 21n|77| as 77 ±00, (9.41) 

and solutions of (9.40) with this asymptotic structure do exist — but not at each end. 
(9.40) admits even solutions, and if we restrict ourselves to these, then we may take 

g'(0) = 0, 3(0)^0. (9.42) 

(If g(0) = 0, then g = is the solution.) However, it is found that such solutions 
have a different asymptotic behaviour as 77 — > 00, namely 



A 

— exp 

M 



(9.43) 



and A = A[g(0)] > for g(0) 7^ (and A(0) = 0), and these cannot match to the 
outer solution. If one alternately prescribes (9.41) as 77 — > +00, for example, then the 
solution is asymmetric, and has the exponential behaviour (9.43) as 77 — > —00. thus 
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the appealingly simple similarity structure implied by steady solutions of (9.39) is 
wrong (and actually, the solution of the initial value problem (9.39) satisfying (9.41) 
tends to zero as r — > oo). 

However, (9.39) itself develops a local similarity structure as r — > oo, using a 
further similarity variable 



z — 



X — Xq 



T l/2 ( to _ t )l/2[_ ln(to _ t )]l/2- 

Rewriting (9.39) in terms of z and r yields 

9t + \zg z + l-e 9 = \\g zz + \zg z \. 
At leading order in r _1 this has a solution 

g = -ln[l + lcz 2 ], 



(9.44) 



(9.45) 



(9.46) 



where c is indeterminate, and this forms the basis for a formal expansion. It is 
algebraically convenient to use (9.46) to define c as a new variable, and also to write 



s = In r; 



Then (9.45) becomes 
2 



TZ* 



2c + 4zc z + z c zz + z 



1 + \cz 2 



C ~~r 2 C 3 



(9.47) 



(9.48) 



We seek a solution for (9.48) in the form 



c ~ cq(z, s) + -ci(z, s) + \c 2 (z, a)..., 



r 



and then, since rd/dr — d/ds, we have 



1/ 1 / 

c s ~ c + -(ci - ci) + — (c 2 - 2c 2 ) + 



(9.49) 



(9.50) 



r t 

where c 2 = dci/ds. Substituting this into (9.48) and equating powers of r, we find 

c = C (s), (9.51) 

where Co is arbitrary, and 



Clz 



2C + z 2 



1 + \W 



+ Co — Co 



(9.52) 



The arbitrary function Co arises because the order of the approximate equation is 
reduced. In order to specify it, and other arbitrary functions of s which arise at each 
order, we require that the solutions q be smooth, and this requires that there be no 
term on the right hand side of (9.52) proportional to 1/z as z — > 0, in order that 
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logarithmic singularities not be introduced. Specifically, we require at each stage of 
the approximation that 



dcj 
~dz~ 



a i + a u z + a 3i z + . . . 



so that z 2 Ci is smooth. Applying this to (9.52) requires that 

Co = Co(l — Co), 
so that Co — > 1 as s — > oo, and then 

2C 



ci = 



+ C 1 ( S ) + C 2 ln[l+|C z 2 ]. 



(9.53) 



(9.54) 



(9.55) 



At 0(1/t 2 ), we then have 



C2z = 



2c 



i + Azc u + z 2 c lzz + z 2 {-(ci - ci) + (ci + |zc ia ) 



2Cq(ci + |ZC i a ) ^ j 



+ \c 2 c 1 z 2 (l + \c z 2 ) 2 



and applying the regularity condition (9.53), we find, after some algebra, 

Ci = 2(1 - C )C! + |C 3 , 



(9.56) 



(9.57) 



so that Ci — >■ Cio + |s as s — > oo. Thus finally we obtain the local similarity solution 



In 



A Ua - t + 



c(x - x ) 2 
4[-\n(t -t)} 



(9.58) 



where c w Co(s), s = lnr = ln[— ln(t — t)]. 

9.4 Notes and references 

There is a large literature on the subject of nonlinear diffusion, blow up, and related 
issues. A thorough examination of some of these problems is contained in the book 
by Samarskii et al. (1995). 



Exercises 

9.1 A small droplet satisfies the surface-tension controlled equation 

h t = -^-V.[/i 3 VV 2 /i]. 

A small quantity / h dV = Q is released at time zero at the origin. Find a 
suitable similarity solution in one and two spatial dimensions. 
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9.2 A gravity-driven drop of fluid spreads out on a flat surface. Its viscosity \x is a 
function of shear rate, so that in (9.1), 

dr 

oz 

(A constant viscosity fluid has n = 1.) Show that the horizontal fluid flux is 

-7U\n— li,n+2i 



|V/i| n_1 /i n+2 V/i, 



n + 2 
and deduce that 

dh = A{pgT v 

dt n + 2 L 1 1 J 
Find similarity solutions in one and two dimensions for the initial emplacement 
of a finite volume at the origin. What happens as n — >• oo or n — > 0? 

9.3 Suppose a two-dimensional drop as in question 9.2 is subjected to a spatially 
varying accumulation a(x), with xd < for x ^ 0. Find appropriate local 
behaviour near the right hand margin x = x s > 0, where h = 0, if x s > 0, 
x s < 0, x s = 0. 

9.4 Let u satisfy 

u t = Xu p + u xx , 

with u = 1 on s = ±1 and t — 0. Prove that if A is large enough, u must blow 
up in finite time if p > 1. Supposing this happens at time to at x — 0, show 
that a possible local similarity structure is of the form 



(tb-t)^ VW ' " (to-*) 1 / 2 ' 
and prove that /3 = (l/(p — 1). Show that in this case, / would satisfy 

f" ~ \W + - Pf = 0, 
and explain why appropriate boundary conditions would be 

/ - l£l" 2/3 as i ±oo, 

and show that such solutions might be possible. Are any other limiting be- 
haviours possible? 

9.5 By direct integration, show that the solution of 

u " + \ e u = 
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satisfying u = Oonx = ±lis 



u = 2 In 




^4sech { \\ 



and find a transcendental equation for A. Hence show that no solution exists 
for A > A c , and derive and solve (numerically) an algebraic equation for A c . 

If the equation is to be solved in [0, 1], with u' = on x = and u' = -1 on 
x — 1, find the solution, and plot u(0) as a function of A. Is there a critical 
value A c ? If so, find it; if not, why not? 

9.6 (i) Find an exact solution of the Gel'fand equation 

V 2 # + \e e = in < r < 1, 

where r is the cylindrical polar radius, and 9 = on r = 1. [Assume cylindrical 
symmetry, and a suitable condition of regularity at r = 0.] Show that there is 
a critical parameter A c such that no solution exists for A > A c , and find its value. 



(ii) Write down the ordinary differential equation satisfied by a spherically sym- 
metric solution of the Gel'fand equation in part (i). Suppose that 9 = on r = 1 
and 9 r = on r = (why?). By putting 

p = \r 2 e e , q = 2 + r9 r , r = e~\ 

show that p(t) and q(t) satisfy the ordinary differential equations 

P = -pq, 

q = p + q — 2. 



By consideration of trajectories for p and q in the phase plane, show that mul- 
tiple solutions exist for A w 2, and infinitely many at A = 2. Sketch the 
corresponding response diagram of 0(0) versus A. 
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